
\part{动态宏观经济学}
\chapter{基本的索罗模型}
经济体只有一种产品，生产函数如，
\[ Y_t=A_t\cdot f(K_t,H_t) \]
其中，$ Y_t $是$ t $期产出，$ A_t $是技术水平，$ K_t $是资本存量，$ H_t $是劳动投入。生产函数满足一次齐次等常见假设，如
\begin{align*}
	F_i(K_t,H_t)>&0\\
	F_{ii}(K_t,H_t)<&0\\
		F_{ij}(K_t,H_t)>&0\\
		\lim_{i\rightarrow 0}F_i(K_t,H_t)\rightarrow& \infty\\
				\lim_{i\rightarrow \infty}F_i(K_t,H_t)\rightarrow& 0\\
				F(K_t,0)=F(0,H_t)=&0, \qquad i=K,H \quad \text{and}\quad i\ne j
\end{align*}

如果$ A_0 $是初期的技术水平，并按$ \alpha $的增长率进行增长，那么$ t $期的技术水平为$ A_t=(1+\alpha)^tA_0 $。同时，为公式中的总量变量进行人均形式的变化，即$ y\equiv Y_t/H_t, k_t\equiv K_t/H_t $，那么依据一次齐次假设，有，
\[ y_t=\frac{Y_t}{H_t}=A_t\cdot f(\frac{K_t}{H_t},\frac{H_t}{H_t})=A_t\cdot f(k_t) \]

设折旧率为$ \delta $，在$ t $期投资为$ I_t $，则资本的增长方程为\footnote{今后我们会明白差分方程是非常重要的。后续的复杂建模中，往往通过最优规划来得到类似的差分方程，从而获得系统的动态演变。但是这里，我们直接给出了差分方程，所以说这个模型是非常简单，非常初步的。}，
\[ K_{t+1}=(1-\delta)K_t+I_t \]
再次通过假设劳动力的增长率\footnote{意味着$ H_t=(1+n)^tH_0 $。}恒为$ n $,那么人均资本的增长方程为，
\[ k_{t+1}=\frac{(1-\delta)k_t+i_t}{1+n}\]
\begin{note}	
注意到上述方程我们引入了一个新的变量——投资，因此，要解这个方程，我们必须增加新的条件。
\end{note}
而这个条件实际上就是模型的均衡条件，即投资要等于储蓄。同时，我们往往令储蓄为产出的一个固定份额，比如$ \sigma $，即
\[ s_t=\sigma y_t \]
那么，人均增长方程可以重新写为，
\begin{equation}\label{solow_eq_a}
(1+n) k_{t+1}=(1-\delta)k_t+\sigma(1+\alpha)^tA_0f(k_t)	
\end{equation}

如果我们认为技术没有进步，即$ \alpha=0 $，那么，上式将更为简单，
\begin{equation}\label{solow_eq_k}
(1+n) k_{t+1}=(1-\delta)k_t+\sigma\cdot A_0f(k_t)	
\end{equation}

上述方程描述了经济的动态演变。同时，有一个术语叫稳态，意味着经济体会稳定不动，此时$ k_{t+1}=k_t=\bar k $，那么\eqref{solow_eq_k}式意味着，
\begin{equation}\label{solow_eq_steady}
(\delta+n)\bar k=\sigma A_0 f(\bar k)	
\end{equation}


这个模型的一个重要结论就是绝对收敛。
\begin{theorem}{绝对收敛}{conv}
不以参数为条件，较贫穷的经济体(拥有较少的资本存量)增长得比富裕经济体(拥有较多初始资本存量)快，最终共同富裕。
\end{theorem}

\begin{proof}
	人均资本增长率为，
	\[\gamma_t=\frac{k_{t+1}}{k_t}=\frac{(1-\delta)k_t+\sigma A_0f(k_t)}{(1+n)k_t}  \]
	该增长率对$ k_t $的导数为，
	\[ \frac{d\gamma_t}{dk_t}=\frac{\sigma A_0}{(1+n)k_t^2}[f'(k_t)k_t-f(k_t)] \]
	注意中括号里的式子在$ k_t>0 $时是小于0的\footnote{中括号里面本质上等价于证明该生产函数的产出弹性小于1，即要素增加1\%，产出增加百分之多少。而一般正常的生产函数都满足产出弹性小于1的性质。回忆高鸿业的《微观经济学》，此处对应于生产的第二阶段，边际产出小于平均产出。
		\begin{align*}
&		f'(k_t)k_t-f(k_t)<0\\
			\Leftrightarrow\;&f'(k_t)<\frac{f(k_t)}{k_t}\\
			\Leftrightarrow\;&\frac{df(k_t)/f(k_t)}{dk_t/k_t}<1
	\end{align*}}。这就意味着，参数一致，贫穷国家增长更快。
\end{proof}
\section{技术进步不为0的情况}
通常地，我们感兴趣稳态。此时的稳态等价于考察内生变量的增长率是一个常数(以前的例子中，稳态意味着增长率为0，但此处是一个常数，意义是一样的)。比如技术进步是一个常数。如果它是一个常数，会发现资本和产出的增长率也是一个常数，而且这个常数与技术进步的增长率是一致的。我们来证明这个结论。

根据\eqref{solow_eq_a}式，如果令$ A_0=1 $，生产函数形式为$ f(k_t)=k_t^\theta $，则增长率可以写为，
\[ \gamma_t=\frac{(1-\delta)k_t+\sigma (1+\alpha)^tk_t^\theta}{(1+n)k_t} \]
上式可以用$ \gamma_t $反算出$ k_t $，
\[ k_t=(1+\alpha)^{\frac{t}{1-\theta}}\left[\frac{\sigma}{(1+n)\gamma_t-(1-\delta)}\right]^{\frac{1}{1-\theta}} \]

因为增长率是常数，所以$ \gamma_t $不用写脚标。上式还意味着，
\[ \gamma=\frac{k_{t+1}}{k_t} =\frac{(1+\alpha)^{\frac{t+1}{1-\theta}}\left[\frac{\sigma}{(1+n)\gamma-(1-\delta)}\right]^{\frac{1}{1-\theta}}}{(1+\alpha)^{\frac{t}{1-\theta}}\left[\frac{\sigma}{(1+n)\gamma-(1-\delta)}\right]^{\frac{1}{1-\theta}}}=(1+\alpha)^{\frac{1}{1-\theta}}\]

而人均产出为，
\begin{align*}
	\frac{y_{t+1}}{y_t}=&\frac{(1+\alpha)^{t+1}k_{t+1}^\theta}{(1+\alpha)^tk_t^\theta}=(1+\alpha)\left(\frac{k_{t+1}}{k_t}\right)^\theta\\
	=&(1+\alpha)\cdot [(1+\alpha)^{\frac{1}{1-\theta}}]^\theta= (1+\alpha)^{\frac{1}{1-\theta}}
\end{align*}

足见，人均资本和人均产出的增长率是相同的，都与技术进步率$ \alpha $相关。
\section{黄金律}
如果消费水平是决定居民福利的唯一因素，那么我们很想知道在稳态时，怎样的储蓄水平最大化了我们的消费？这个储蓄水平就是黄金储蓄规则(Golden rule)。

稳态时的消费可以表达为，
\begin{equation}\label{solow_eq_c}
\bar{c}=(1-\sigma)\bar{y}=(1-\sigma)A_0f(\bar k)	
\end{equation}


而$ \sigma $与$ \bar k $是有函数关系的，不是孤立的，所以上式针对$ \sigma $(或者$ \bar k $)最大化时，要把这个函数关系代进去。这个函数关系就是前面关于稳态的\eqref{solow_eq_steady}式，将其带入\eqref{solow_eq_c}式，有，
\[ \bar c = A_0f(\bar k)-(\delta+n)\bar k \]
上式关于$ \bar k $最优化的一阶条件表明，
\[ A_0f'(\bar k^*)=\delta+n \]

上式可以得到最优的$ \bar k $，即$ \bar k^* $,然后用$ \bar k^* $代入稳态方程就可以得到黄金储蓄率，
\[ \sigma_{GR}=\frac{(\delta+n)\bar k^*}{A_0f(\bar k^*)} \]

\section{随机Solow模型}
现在将模型进一步接近实际，加入随机冲击，比如令某变量符合某个简单的随机过程。这个变量可以是技术水平，贴现因子，储蓄率，人口增长率等。通常我们会对技术水平随机化。

如果令技术水平服从一个AR(1)过程，如
\[ A_t=\psi \bar A +(1-\psi)A_{t-1}+\varepsilon_t \]
其中$ \bar A $是常数，$ \psi $是系数，则随机冲击项$ \varepsilon_t $必须要大于$ -\psi \bar A $。否则技术水平可能为负。为了避免这个麻烦，我们可以用下面这种随机过程，
\begin{equation}\label{solow_eq_At}
A_t=\bar A e^{\varepsilon_t}	
\end{equation}

这就意味着，
\[ \ln A_t=\ln \bar A+\varepsilon_t \]

而索罗模型的一阶差分方程为\eqref{solow_eq_k}式，该式在技术增长遵循\eqref{solow_eq_At}式时的形式为，
\begin{equation}\label{solow_eq_sto_k}
k_{t+1}=\frac{(1-\delta)k_t+\sigma\cdot\bar Ae^{\varepsilon_t}f(k_t)}{1+n}	
\end{equation}


\begin{note}
	在动态建模中，一阶差分方程就是我们模型的解。有了这个方程，我们才可以算稳态，算后续的各种脉冲响应。
	\end{note}

利用\eqref{solow_eq_sto_k}是可以模拟经济中人均资本的动态变化情况的。比如在$ \bar A=1, n = 0.02, \delta=0.1,\theta=0.36,\sigma=0.2,\varepsilon_t $的标准差为0.2时，我们随机模拟了三次如下图所示，
\begin{figure}[H]
	\includegraphics[scale=0.7]{solow1.png}
	\caption{三次模拟}
	源码文件\lstinline|simu_solow.R|.
\end{figure}

\begin{lstlisting}[language=R,mathescape=false,numbers=left]
rm(list = ls())
library(tidyverse)
library(ggplot2)
# parameters set
A <- 1
dlt <- 0.1
sgm <- 0.2
theta <- 0.36
n <- 0.02
eps <- rnorm(120, sd = 0.2)

# 迭代. 公式(1.2)
compk <- function(A, dlt, sgm, theta, n, eps){
	k <- rep(2.2,120) # 初值2.2
	for (i in 2:120) {
		k[i] <- ((1-dlt)*k[i-1] + sgm*A*exp(eps[i-1])*k[i-1]^theta)/(1+n)
	}
	return(data.frame(k = k))
}

# 模拟三次
set.seed(1024)
eps <- lapply(rep(0.2,3),rnorm, n = 120, mean = 0)
picdata <- lapply(eps, compk, A = A, dlt = dlt, sgm = sgm, theta = theta, n = n) %>% bind_cols()
names(picdata) <- c('k1','k2','k3')
picdata$time <- 1:120

# 画图
ggplot(picdata, aes(x = time, y = k1)) + geom_line(linetype = 1) +
geom_line(aes(y = k2), linetype = 2) +
geom_line(aes(y = k3), linetype = 9) + labs(y = 'k') +
theme_bw()
\end{lstlisting}

\section{对数线性形式的Solow模型}
但问题在于，你可以模拟$ k_t $，也可以模拟$ y_t $，而实际上$ y_t $是可以观察的，你怎么知道你用上述参数模拟的$ y_t $的方差就符合现实中的数据呢？

为了从冲击的方差可以解析出$ k_t $或$ y_t $的方差，对于我们前面的非线性模型是不可能的，必须在稳态附近对模型进行线性近似。

为了对数线性化，首先把非线性方程写成对数差分形式。具体如下，定义
\[\widetilde X_t=lnX_t-ln\overline X\]

就有\footnote{我们对非线性方程通常进行这种对数差分变化。}，
\begin{equation}\label{sl_eq_x}
X_t=\overline Xe^{\widetilde X_t}	
\end{equation}

然后就可以利用下面这些等式对\eqref{sl_eq_x}式进行线性化\footnote{第一个等号成立的源自在于指数函数在原点的泰勒展开。即，\[ e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots \]此处仅利用了一阶近似。}，
\begin{align*}
	e^{\widetilde X_t}&\approx 1+\widetilde X_t\\
	\widetilde X_t\widetilde Y_t&\approx 0\hspace{2em}\text{高阶的消去而已}\\
	E_t[ae^{\widetilde X_{t+1}}]&\approx a+aE_t(\widetilde X_{t+1})	\hspace{2em}\text{第一个法则在期望里面展开而已}\\
\end{align*}

\subsection{资本的方差}
我们的目标在于获得$ k_t $与$ \varepsilon_t $之间的函数关系。

对于一阶差分方程$(1+n) k_{t+1}=(1-\delta)k_t+\sigma\cdot\bar Ae^{\varepsilon_t}f(k_t) $，将$ k_j $换成$ \bar ke^{\widetilde{k}_j}  $，得到\footnote{我们利用了$ f(k_t)=k_t^\theta $。}，
\[ (1+n)\bar ke^{\widetilde{k}_{t+1}}=(1-\delta)\bar ke^{\widetilde{k}_t}+\sigma\cdot\bar Ae^{\varepsilon_t}\bar k^\theta e^{\theta\widetilde{k}_t}	 \]
再利用$ e^{\widetilde X_t}\approx 1+\widetilde X_t $，则上式可以写为，
\[ (1+n)\bar k(1+\widetilde{k}_{t+1})=(1-\delta)\bar k(1+\widetilde{k}_{t})+\sigma\cdot\bar A(1+\varepsilon_t)\bar k^\theta (1+\theta\widetilde{k}_{t})	 \]

化简上面这个式子，就能得到$ \widetilde{k}_t$与$ \varepsilon_t $的函数关系。在化简时，如果用到稳态时的关系\eqref{solow_eq_steady}式会有助于化简，同时也用到$ \widetilde X_t\widetilde Y_t\approx 0 $此处体现为$ \varepsilon_t\tilde{k}_t\approx 0 $。有，
\begin{equation}\label{solow_eq_rec_k}
\widetilde{k}_{t+1}=B\widetilde{k}_t+C\varepsilon_t 	
\end{equation}

其中，
\begin{align*}
	B = & \frac{1+\theta n-\delta(1-\theta)}{1+n}<1\\
	C = &\frac{\delta+n}{1+n}
\end{align*}
对于\eqref{solow_eq_rec_k}式，可以递归得到，
\[ \widetilde{k}_{t+1}=C\sum_{i=0}^\infty B^i\varepsilon_{t-i} \]
基于上式，计算$ \widetilde{k}_{t+1} $的方差，有，
\[ var(\widetilde{k}_{t+1}) = E(\widetilde{k}_{t+1}\widetilde{k}_{t+1}) = E\left[\left(C\sum_{0}^\infty B^i\varepsilon_{t-i}\right)\left(C\sum_{0}^\infty B^i\varepsilon_{t-i}\right)\right]\]
如果冲击独立，则上式等于，
\[ var(\widetilde{k}_{t+1})=C^2\sum_{i=0}^\infty B^{2i}var(\varepsilon)=\frac{C^2}{1-B^2}var(\varepsilon),\;\;\text{等比数列求和} \]
其中，
\[ \frac{C^2}{1-B^2}=\frac{(\delta+n)^2}{(1+n)^2+[1+\theta n-\delta(1-\theta)]^2} \]
利用前面参数的取值，可知资本存量的方差是技术冲击方程的0.0955倍。

\begin{figure}[H]
\includegraphics[scale=0.7]{solow2.png}
\caption{解析解和对数线性化解的模拟对比}
源码文件\lstinline|simu_solow.R|.
\end{figure}
\begin{lstlisting}[language=R,numbers=left,firstnumber=33,mathescape=false]
# 对数线性化的迭代. 公式(1.4)
B <- (1+theta*n-dlt*(1-theta))/(1+n)
C <- (dlt+n)/(1+n)

compk_ln <- function(B, C, eps){
	k <- rep(0,120) # 因为是对数差分的情况，所以初值为0
	for (i in 2:120) {
		k[i] <- B*k[i-1] + C*eps[i-1]
	}
	return(data.frame(k = k))
}


# 画图
picdata$k_ln <- compk_ln(eps = eps[[1]],B = B, C = C) %>%
unlist()
picdata$k_ln <- 2.2*exp(picdata$k_ln) # 由对数差分返回原值

ggplot(picdata, aes(x = time, y = k1)) + geom_line(linetype = 1) +
geom_line(aes(y = k_ln), linetype = 2) + labs( y = 'k') +
theme_bw()
\end{lstlisting}

\subsection{产出的方差}
\begin{align*}
	y_t=&\bar A e^{\varepsilon_t}k_t^\theta\\
	\Rightarrow \qquad \bar ye^{\widetilde{y}_t}=&\bar Ae^{\varepsilon_t}\bar k_t^\theta e^{\theta\widetilde{k}_t}\qquad\text{稳态替换}\\
	\Rightarrow\qquad\bar y(1+\widetilde{y}_t)=&\bar A(1+\varepsilon_t)\bar k^\theta (1+\theta \widetilde{k}_t)\qquad\text{近似替换}
\end{align*}
化简，得，
\[ \widetilde{y}_t=\varepsilon_t+\theta \widetilde{k}_t \]

因为$ \widetilde{k}_t $独立于$ \varepsilon_t $，故，
\[ var(\widetilde{y}_t)=E(\widetilde{y}_t\widetilde{y}_t)=E[(\varepsilon_t+\theta \widetilde{k}_t)(\varepsilon_t+\theta \widetilde{k}_t)]=var(\varepsilon_t)+\theta^2var(\widetilde{k}_t) \]

可以得到：
\begin{itemize}
	\item $ var(\widetilde{y}_t)=1.0123var(\varepsilon_t) $,产出方差与技术冲击方差几乎一致；
	\item $ var(\widetilde{y}_{t+1})=var(\varepsilon_{t+1}) +\theta^2\frac{C^2}{1-B^2}var(\varepsilon_t)$，而$ \theta^2\frac{C^2}{1-B^2}=0.0123 $，意味着$ t $期的冲击只能解释$ t+1 $期产出波动的1.23\%，技术冲击的持久性太弱。
\end{itemize}
\begin{remark}
上述两点意味着，我们的参数设置可能是不合理的。
\end{remark}

\chapter{OLG模型的储蓄}
索罗模型的一个大的缺陷就是储蓄率是外生决定的，现在我们用世代交叠模型来内生化储蓄率(当然还有其他内生化的方法，下一章讲述)。

\paragraph{故事的梗概：}
\begin{itemize}
	\item 模型中的个体生命期为2期，第一期是年轻时期，第二期是退休的年老时期；
	\item 行为人应在第一期对其两期的生命消费做出规划；
	\item 每个人出生在当期，然后工作、消费直到第2期老去。这个过程无限延续。
\end{itemize}


\paragraph{符号约定：}
\begin{table}[H]
\begin{tabular}{lp{20em}}\hline
符号&含义\\\hline
$ N(t) $&第$ t $代有$ N(t) $个个体\\
$ u_t(c_t^h(t),c_t^h(t+1)) $&第$ t $代第$ h $个个体的效用函数。脚标$ t $表示出生的时期。要注意脚标$ t $和括号里面的$ t $的不同含义。$ c_t^h(s) $指的是第$ t $代第$ h $个人在时期$ s $的消费。\\
$Y(t)=F(K(t),H(t))$&生产函数\\
$ H(t) $&劳动总量\\
$ h_t^h=[h_t^h(t),h_t^h(t+1)] $&个体劳动禀赋由2期构成\\
$ H(t)=\sum_{h=1}^{N(t)}h_t^h(t)+\sum_{h=1}^{N(t-1)}h_{t-1}^h(t) $&劳动总量构成\\
$ H_t(t)= \sum_{h=1}^{N(t)}h_t^h(t)$&$ t $期年轻人总劳动禀赋\\
$ H_{t-1}(t)= \sum_{h=1}^{N(t-1)}h_{t-1}^h(t)$&$ t $期老年人总劳动禀赋\\\hline
\end{tabular}
\end{table}



每个个体的最优规划为，
\begin{align*}
	\max_{c_t^h(t),c_t^h(t+1)}\qquad&u_t^h[c_t^h(t),c_t^h(t+1)]\\
	s.t.\qquad& c_t^h(t)+\frac{c_t^h(t+1)}{r_{t+1}}=w_th_t^h(t)+\frac{w_{t+1}h_t^h(t+1)}{r_{t+1}}
\end{align*}

该约束条件的经济含义为一生总消费的现值等于其收入的现值。可以取个巧，将预算约束代入上述最优规划，使得上述规划变成无约束规划，
\[ \max_{c_t^h(t)}\qquad u[c_t^h(t), r_{t+1}w_th_t^h(t)+w_{t+1}h_t^h(t+1)-r_{t+1}c_{t}^h(t)] \]
上述最优规划的一阶条件为，
\begin{equation}\label{olg_eq_u}
u_1=r_{t+1}u_2	
\end{equation}

其中，$ u_i $表示对$ u $的第几个变量求偏导。注意，上式隐含地表达了$ c_t^h(t) $是$ r_t,w_t,w_{t+1} $的函数(劳动禀赋我们通常将其当做一个参数看待)。同时，我们注意到系统还有一个约束，即在$ t $期，第$ t $代人工作而获得收入，该收入用来在$ t $期消费，并进行下一期的资本累积，即年轻时有预算约束，
\begin{equation}\label{olg_eq_young}
	w_th_t^h(t)=c_t^h(t)+k^h(t+1)
\end{equation}


将\eqref{olg_eq_u}式隐含的$ c_t^h(t) $带入上式，然后累加，就会得到，
\begin{equation}\label{olg_eq_bug}
 w_tH_t(t)=\sum_{h=1}^{N(t)}c_t^h(t)+K(t+1) 	
\end{equation}


同时注意到$ w_t=F_H[K(t),H(t)],r_t= F_K[K(t),H(t)]$，你会发现\eqref{olg_eq_bug}式是一个关于$ K(t),K(t+1) $的差分方程，这就是我们要的解。

另外，可以看一下上述最优规划的约束条件如何得到。因为年老时的预算约束为，
\begin{equation}\label{olg_eq_old}
c_t^h(t+1)=w_{t+1}h_t^h(t+1)+r_{t+1}k^h(t+1)	
\end{equation}

那么结合\eqref{olg_eq_young}式消去$ k^h(t+1) $，就会有人一生的预算约束为，
 \[ c_t^h(t)+\frac{c_t^h(t+1)}{r_{t+1}}=w_th_t^h(t)+\frac{w_{t+1}h_t^h(t+1)}{r_{t+1}} \] 
 
 
\section{一个实例}
假设消费者最优规划为，
\begin{align*}
	\max_{c_t^h(t),c_t^h(t+1)}\qquad&u_t^h[c_t^h(t),c_t^h(t+1)]=c_t^h(t)c_t^h(t+1)^\beta\\
	s.t.\qquad& c_t^h(t)+\frac{c_t^h(t+1)}{r_{t+1}}=w_th_t^h(t)+\frac{w_{t+1}h_t^h(t+1)}{r_{t+1}}
\end{align*}
则其一阶条件对应为，
\begin{align*}
c_t^h(t+1)^\beta = &r_{t+1}\beta c_t^h(t)c_t^h(t+1)^{\beta-1} \\
\Rightarrow\qquad c_t^h(t+1) = &r_{t+1}\beta c_t^h(t)
\end{align*}
然后利用年轻和年老的约束，即\eqref{olg_eq_young}式和\eqref{olg_eq_old}式，汇总到总量层面，有，
\begin{equation}\label{olg_eq_K1}
  w_{t+1}H_t(t+1)+r_{t+1}K(t+1)=r_{t+1}\beta [w_tH_t(t)-K(t+1)]
\end{equation}

如果生产函数形式为$ Y_t=K(t)^\theta H(t)^{1-\theta} $，则工资和利率为，
\[ w_t= (1-\theta)K(t)^\theta H(t)^{-\theta},\qquad r_t= \theta K(t)^{\theta-1}H(t)^{1-\theta}\]
将工资和利率带入\eqref{olg_eq_K1}式，整理，有，
\begin{align*}
K_{t+1}=\frac{\theta\beta\frac{H_t(t)}{H(t)^\theta}}{\frac{H_t(t+1)}{H(t+1)}+\frac{\theta(1+\beta)}{1-\theta}}\cdot K(t)^\theta=\kappa \cdot K(t)^\theta
\end{align*}
可以参见Maple软件的符号运算\lstinline|solve_K.mw|。这个模型的关键就是建立最优规划以后，利用约束将$ K(t) $和$ K(t+1) $联系起来。

如果取$ \kappa =4,\theta = 0.36 $，我们也可以模拟该方程的动态演化如下所示，稳态在曲线与45\textdegree 线的交点处。


如果定义总资本增长率为，
\[ \Delta K(t)=\frac{K(t+1)}{K(t)}= \frac{\kappa K(t)^\theta}{K(t)}\]
那么它对资本的导数为，
\[ \frac{d\Delta K(t)}{dK(t)}=(\theta-1)\kappa K(t)^{\theta-2}<0 \]
即资本存量越大，增长速度越慢，依然是绝对收敛。类似地，我们也可以导出产出增长率对资本的导数，
\[ \Delta Y(t)=\frac{K(t+1)^\theta H(t+1)^{1-\theta}}{K(t)^\theta H(t)^{1-\theta}}=\Delta K(t)^\theta \]

因为劳动供给是常数，所以$ H(t)=H(t+1) $。产出增长率对资本存量的导数为，
\[ \frac{d\Delta Y(t)}{dK(t)}=\frac{d\left[\frac{\kappa K(t)^\theta}{K(t)}\right]^\theta}{dK(t)}=\theta\left[\frac{\kappa K(t)^\theta}{K(t)}\right]^{\theta-1}(\theta-1)\kappa K(t)^{\theta-2}<0 \]
\section{随机版本}
生产函数为，
\begin{align*}
Y_t=&\lambda_tK(t)^\theta H(t)^{1-\theta}\\
\lambda_t=&(1-\gamma)+\gamma \lambda_{t-1} +\varepsilon_t,	\qquad0<\gamma<1,E_t(\varepsilon_{t+1})=0
\end{align*}
厂商利润最大化要求，
\begin{align*}
	w_t=&(1-\theta)\lambda_tK(t)^\theta H(t)^{-\theta}\\
	r_{t+1}=&\theta\lambda_{t+1}K(t+1)^{\theta-1}H(t+1)^{1-\theta}
\end{align*}

消费者最优化问题为(注意到这里的效用函数与前面的不一样，这里取了对数)，
\begin{align}
	\max_{c_t^h(t),c_t^h(t+1)} \quad & E_t[u(c_t^h(t),c_t^h(t+1))]=\ln c_t^h(t)+\beta E_t\ln c_t^h(t+1)\\\label{olg_eq_y2}
	s.t.\qquad & w_th_t^h(t)=c_t^h(t)+k^h(t+1)\qquad\text{年轻人约束}\\\label{olg_eq_old2}
	&c_t^h(t+1)=w_{t+1}h_t^h(t+1)+r_{t+1}k^h(t+1)\qquad\text{老人约束}	
\end{align}	
	
	可以联合两个约束条件，替换掉目标函数中的$ c_t^h(t+1) $，然后变成单变量优化问题。即，
	\[ \max_{c_t^h(t)}\quad\ln c_t^h(t)+\beta E_t\ln [r_tw_th_t^h(t)+w_{t+1}h_t^h(t+1)-r_{t+1}c_t^h(t)] \]
	
	该问题的一阶条件为，
	\[ \frac{1}{c_t^h(t)}=\beta E_t\left[\frac{r_{t+1}}{c_{t}^h(t+1)}\right] \]
	
	将年轻人约束\eqref{olg_eq_y2}式和老人约束\eqref{olg_eq_old2}式汇总后，代入上式，替换掉$ c_t^h(t),c_t^h(t+1) $，可以得到，
	\begin{equation}\label{olg_eq_K}
K(t+1)=\kappa \lambda_tK(t)^\theta		
	\end{equation}
	其中，
	\[ \kappa = \frac{\theta\beta\frac{H_t(t)}{H(t)^\theta}}{\frac{H_t(t+1)}{H(t+1)}+\frac{\theta(1+\beta)}{1-\theta}} \]
	与前面的模型是一样的。
	
当给定参数值$ \beta=0.99,\theta=0.36,H(t)=65,H_t(t)=60.2,H_t(t+1)=4.8 $时，$ \kappa=4 $，再假设$ \gamma=0.6 $，$ \varepsilon_t $在[-0.02,0.02]上服从均匀分布，那么以$ K(0)=(6.9792,8.7241,10.4689) $为初始值，根据\eqref{olg_eq_K}式，就可以模拟$ K $的走势，根据生产函数就可以模拟出产出走势，走势如下图所示，

\begin{figure}[H]
\includegraphics[scale=0.7]{olg.png}
\caption{OLG模型的三次不同初始值下的模拟}
\small 本图源码为\lstinline|simu_olg.R|
\end{figure}


%\paragraph{故事的解释：}
%$ t $期的资本存量$ K(t) $是从$ t-1 $期继承过来的，而且在$ t $期不能改变，而且这部分资本会在$ t $期完全使用掉。这就意味着：
%\[ Y(t)=F[K(t),H(t)]\ge \sum_{h=1}^{N(t)}c_t^h(t) +\sum_{h=1}^{N(t-1)}c_{t-1}^h(t)+K(t+1)\]
%即本期的产出为年轻人和老人所消费，以及留给下一期的资本存量。

\chapter{无限期生存的行为人}\label{inf}
另一种消除储蓄率外生的办法。变分法求得稳态：根本思路在于先假定内生变量在$ s-1 $和$ s+1 $期的值都是给定的，再利用一阶条件最大化$ s $期的目标方程，然后因为稳态时内生变量各期值应该相同，就此就可以得到稳态。
\section{固定劳动投入的RC经济}
举例最为方便。在RC模型中，时期$ t $的消费为，
\[ c_t=f(k_t)-k_{t+1}+(1-\delta)k_t \]

它是由资本积累方程$ k_{t+1}=(1-\delta)k_t+i_t$和$ y_t=c_t+i_t $而得到的。毕生效用函数为，
\[ \sum_{t=0}^\infty \beta^t u[f(k_t)-k_{t+1}+(1-\delta)k_t] \]

\textbf{把$ k_t,k_{t+2} $均看作给定的，则可以得到关于$ k_{t+1} $期的一阶条件为}\footnote{本质上，这里要用到我们前面讲过的动态优化理论，但你会发现你构造汉密尔顿函数，然后巴拉巴拉，会得到我们这种简单处理的一样的结果。因此，以后像这样的动态优化，都可以这样简单快捷地处理。}，
\[ u'[f(k_t)-k_{t+1}+(1-\delta)k_t]=\beta u'[f(k_{t+1})-k_{t+2}+(1-\delta)k_{t+1}][f'(k_{t+1})+(1-\delta)] \]

该一阶条件我在\lstinline|RC_oneorder.mw|中有一个简单的计算。再令$ k_{t+1}=k_{t+2}=k_t=\bar k $，则可以得到，
\[ f'(\bar k)=\frac{1}{\beta}-1+\delta \]

该一阶条件只是必要条件，为了确定最终的解，通常还有一个横截条件或者说边界条件。边界条件的一般形式是，
\[ \lim_{t\rightarrow\infty}\beta^tu_1(k_t,k_{t+1})k_t=0 \]
含义为在足够远的未来，资本存量增加一个单位带来的边际效用贴现到现在的值为0。否则，无限期地推迟消费将是最优解，稳态也不会存在。可以考察一个横截条件不成立的例子。

\paragraph{横截条件不成立的例子}考虑一个经济体，$ f(k_t)=0.2k_t,u(c_t)=c_t,\delta=0.1,\beta=0.98 $。假设个体将在$ T $期死亡。其最大化目标函数为，
\begin{align*}
	\max\;&\sum_{t=0}^T\beta^tu[f(k_t)-k_{t+1}+(1-\delta)k_t] =\sum_{t=0}^{T}0.98^t(1.1k_t-k_{t+1})\\
	s.t.\;\;&c_t=1.1k_t-k_{t+1}\ge 0
\end{align*}
拉格朗日函数为，
\[ L= \sum_{t=0}^{T}[0.98^t(1.1k_t-k_{t+1})+\lambda_t(1.1k_t-k_{t+1})]\]
基于库恩塔克条件，一阶条件如下，
\begin{align*}
	\frac{\partial L}{\partial k_{t+1}}=-0.98^t+0.98^{t+1}\cdot 1.1-\lambda_t+1.1\lambda_{t+1}=0\\
	c_t=1.1k_t-k_{t+1}\ge 0,\;\; \lambda_t\ge 0\text{互补松弛}
\end{align*}
根据第一个式子可知$ \lambda_t,\lambda_{t+1} $不可能为0，故互补松弛条件中关于$ c_t $的约束条件的等号必然成立，即每期的消费为0，每期的资本存量都是上一期资本存量的1.1倍，直到最后一期一次性消费全部累积的资产（注意$ k_{T+1}=0 $）。若给定初值$ k_0 $，此时最大化的效用就是角点解，
\[ (0.98\times 1.1)^Tk_0 \]
若在无限期情况下，效用最大化的途径是无限推迟消费，目标函数最大值是无穷远处资本存量效用贴现值，因为此时的横截条件是，
\[ \lim_{t\rightarrow\infty}(0.98^t\times 1.1)\cdot 1.1^{t-1}k_0 \]
而它是不等于0的。我们要避免出现这种情况。尽管我们不是每次都检验这个条件，但研究者一定要确保横截条件成立。一般地，当效用函数和生产函数都是凹的、连续可微的、各期的可行集都是凸的和非空的，那么该经济体是标准的，横截条件是成立的。

\section{可变劳动投入的RC经济}
为了更加接近现实，可以允许行为人决定如何在闲暇和工作之间分配时间。此时，个体效用函数为，
\begin{align}\nonumber
\max_{c_t,h_t}\quad\sum_{t=0}^\infty & \beta^t\bar u(c_{t},l_{t})=\sum_{t=0}^\infty \beta^t \bar{u}(c_{t},1-h_{t})=\sum_{t=0}^\infty \beta^t u(c_{t},h_{t})\\\label{vl1}
s.t.\hspace{2em} & k_{t+1}=(1-\delta)k_t+i_t\\\label{vl2}
& y_t=f(k_t,h_t)=c_t+i_t
\end{align}
其中，$ l_t $是闲暇，从而$ h_t=1-l_t $是工作时间。人们需要选择消费$ c_t $、工作时间$ h_t $，投资$ i_t $以及资本$ k_{t+1} $来获得最大的效用。

该最优规划的拉格朗日函数为\footnote{也可以把约束条件都代入目标函数中，使得仅对$ k_{t+1},h_t $进行优化。如
\[ \max_{k_{t+1},h_t}\quad \sum_{t=0}^\infty \beta^tu[f(k_t,h_t)-k_{t+1}+(1-\delta)k_t,h_t] \]
其一阶条件可以迅速得到，
\begin{align*}
	\frac{u_h(c_t,h_t)}{u_c(c_t,h_t)}&=-f_h(k_t,h_t)\\
	\frac{u_c(c_t,h_t)}{u_c(c_{t+1},h_{t+1})}&=\beta[f_k(k_{t+1},h_{t+1})+(1-\delta)]
\end{align*}

}，
\[ \mathcal{L}=\sum_{i=0}^\infty \beta^t[u(c_{t},h_{t})-\lambda_{1,t}(k_{t+1}-(1-\delta)k_t-i_t)-\lambda_{2,t}(f(k_t,h_t)-c_t-i_t)] \]
四个一阶条件为，
\begin{align}\label{vl3}
\frac{\partial \mathcal{L}}{\partial c_t}& = u_c(c_t,h_t)+\lambda_{2,t}=0\\\label{vl4}
\frac{\partial \mathcal{L}}{\partial h_t}& = u_h(c_t,h_t)-\lambda_{2,t}f_h(k_t,h_t)=0\\\label{vl5}
\frac{\partial \mathcal{L}}{\partial k_{t+1}}& =-\lambda_{1,t}+\beta\lambda_{1,t+1}(1-\delta)-\beta\lambda_{2,t+1}f_k(k_{t+1},h_{t+1})=0\\\label{vl6}
\frac{\partial \mathcal{L}}{\partial i_{t}}&=\lambda_{1,t}+\lambda_{2,t}=0
\end{align}

方程\eqref{vl1}式至\eqref{vl6}式构成整个系统。先整理一下，以消去两个$ \lambda $，
\begin{align}\label{vl7}
\frac{u_h(c_t,h_t)}{u_c(c_t,h_t)}&=-f_h(k_t,h_t)\\\label{vl8}
\frac{u_c(c_t,h_t)}{u_c(c_{t+1},h_{t+1})}&=\beta[f_k(k_{t+1},h_{t+1})+(1-\delta)]
\end{align}

此时可以考察均衡状态，即忽略时间脚标，那么首先是资本存量的\eqref{vl8}式可以写成，
\begin{equation}\label{vl9}
 \frac{1}{\beta}-(1-\delta)=f_k(k,h)
\end{equation}
\eqref{vl7}式写成，
\begin{equation}\label{vl10}
 \frac{u_h(c,h)}{u_c(c,h)}=-f_h(k,h)
\end{equation}
\eqref{vl1}式写成，
\begin{equation}\label{vl11}
 \delta k=f(k,h)-c
\end{equation}

三个未知数$ c,h,k $，三个方程，\eqref{vl9}-\eqref{vl11}。
\subsection{一个实例}
如果令生产函数为
\[ y_t=k_t^\theta h_t^{1-\theta} \]
效用函数为
\[ u(c_t,h_t)=\ln c_t+B\ln (1-h_t) \]
那么，\eqref{vl9}-\eqref{vl11}式对应的方程为，
\begin{align*}
	 \frac{1}{\beta}-(1-\delta)=&\theta\bar{k}^{\theta-1}\bar{h}^{1-\theta}\\
	 \frac{B\bar c}{1-\bar{h}}=&(1-\theta)\bar{k}^\theta \bar{h}^{-\theta}\\
	 \delta \bar{k}=&\bar{k}^\theta\bar{h}^{1-\theta}-\bar c
\end{align*}
三个方程，三个未知数，我们很容易得到解，文件\lstinline|RC.mw|给出了解析解。在真实的求解中，我们很少会得到解析解，我们一般会数值求解非线性方程组。

\section{竞争经济}
前面考察的是只有一个人的经济，现在假设许多完全相同的个体构成行为人总体，并标准化为1，即总人数是1。消费者向市场供给劳动，厂商以竞争性工资$ w_t $雇佣劳动，其生产函数为$ f(K_t,H_t) $，大写表示社会中的资本和劳动总量，其使用资本的利率为$ r_t $。因为厂商是竞争的，所以劳动边际产出等于工资，资本边际产出等于利率。

个体$ i $的最优规划为，
\begin{align*}
\max&\sum_{t=0}^\infty \beta^tu(c_t^i,h_t^i)\\
s.t.\hspace{1em}&c_t^i=w_th_t^i+r_tk_t^i-I_t^i\\
& w_t=f_h(K_t,H_t)\\
& r_t=f_k(K_t,H_t)\\
&k_{t+1}^i=(1-\delta)k_t^i+I_t^i
\end{align*}
$ I_t^i $表示个体$ i $在$ t $期的投资。那么把工资和利息方程代入消费约束，则上式的拉格朗日函数可以写成，
\[ \mathcal{L}^i=\sum_{t=0}^\infty\beta^t[u(c_t^i,h_t^i)-\lambda_{1,t}(k_{t+1}^i-(1-\delta)k_t^i-I_t^i)-\lambda_{2,t}(f_h(K_t,H_t)h_t^i+f_k(K_t,H_t)k_t^i-c_t^i-I_t^i)] \]

类似的，四个一阶条件为，
\begin{align*}
\frac{\partial \mathcal{L}^i}{\partial c_t^i}& = u_c(c_t^i,h_t^i)+\lambda_{2,t}=0\\
\frac{\partial \mathcal{L}}{\partial h_t^i}& = u_h(c_t^i,h_t^i)-\lambda_{2,t}f_h(K_t,H_t)=0\\
\frac{\partial \mathcal{L}}{\partial k^i_{t+1}}& =-\lambda_{1,t}+\beta\lambda_{1,t+1}(1-\delta)-\beta\lambda_{2,t+1}f_k(K_{t+1},H_{t+1})=0\\
\frac{\partial \mathcal{L}}{\partial I^i_{t}}&=\lambda_{1,t}+\lambda_{2,t}=0
\end{align*}
消去两个$ \lambda $，再联合预算约束，得到，
\begin{align}\label{vl13}
\frac{u_h(c_t,h_t)}{u_c(c_t,h_t)}&=-f_h(K_t,H_t)\\\label{vl14}
\frac{u_c(c_t,h_t)}{u_c(c_{t+1},h_{t+1})}&=\beta[f_k(K_{t+1},H_{t+1})+(1-\delta)]\\\label{vl12}
k_{t+1}^i&=(1-\delta)k_t^i+f_h(K_t,H_t)h_t^i+f_k(K_t,H_t)k_t^i-c_t^i
\end{align}
以及两个加总规则，
\[ H_t=\int_{0}^{1}h_t^idi,\hspace{2em}K_t=\int_{0}^1k_t^idi \]
注意这个积分，如果个体都是一模一样的，那么这个积分里面的这个函数，$ h_t^i $或者$ k_t^i $就是个常数，对这个常数在[0,1]区间定积分，很明显，
\[ H_t=h_t^i,\hspace{2em}K_t=k_t^i \]

如果生产函数一次齐次(规模报酬不变)，厂商完全竞争，从而零利润。这就意味着下式成立\footnote{对于厂商的最优规划可以写为，\[ \max_{K_t,H_t}f(K_t,H_t)-rK_t-wH_t \]很明显其一阶条件为，\begin{align*}
		f_{K_t}=&r\\ f_{H_t}=&w
\end{align*}
而且竞争均衡下，利润为0，那么将一阶条件带入目标函数，有\[ f_k(K_t,H_t)H_t+f_k(K_t,H_t)K_t=f(K_t,H_t)  \]}，
\[ f_k(K_t,H_t)H_t+f_k(K_t,H_t)K_t=f(K_t,H_t) \]
如果我们把预算约束\eqref{vl12}式进行个体汇总，再利用上式就有，
\[ K_{t+1}=(1-\delta)K_t+f(K_t,H_t)-\int_{0}^1c_t^idi=(1-\delta)K_t+f(K_t,H_t)-C_t \]

因为个体相同，一阶条件可以直接替换成总量形式，
\begin{align}\label{vl15}
\frac{u_h(C_t,H_t)}{u_c(C_t,H_t)}&=-f_h(K_t,H_t)\\\label{vl16}
\frac{u_c(C_t,H_t)}{u_c(C_{t+1},H_{t+1})}&=\beta[f_k(K_{t+1},H_{t+1})+(1-\delta)]\\\label{vl17}
K_{t+1}&=(1-\delta)K_t+f(K_t,H_t)-C_t
\end{align}

比较个体均衡的\eqref{vl13}式至\eqref{vl12}式，以及总体均衡的\eqref{vl15}式至\eqref{vl17}，会发现它们竟然是一模一样的。

要注意观察RC经济(也可以称为中央计划者经济)的最优优化与竞争经济的最优优化，前者并不考虑价格(工资和利率)。

\section{第二福利定理}

\begin{theorem}{福利经济学第二定理}
	如果消费者的偏好和厂商的生产集是凸的，则市场是完全的，且市场价格是众所周知的，每个行为人都是价格接受者，那么，只要允许安排适当的一次性财富转移支付，则任何帕累托最优结果都可以通过竞争均衡来实现。
\end{theorem}

这里举一个包含扭曲的例子，说明此时的竞争经济与上一小节讲的代表性行为人经济的均衡是不一样的。这也意味着如果存在扭曲，就不能用代表性行为人去模仿这个经济。在模型中加入政府，政府征收工资所得税，然后通过一次性转移支付把税再返还给每个家庭。\textbf{这样的设置，不会使得代表行为人的均衡条件出现任何变化，但竞争性经济的均衡条件则会发生改变。}

设税率为$ t_w $，$ t $期一次性转移的总量为$ T_t $，因此，政府的预算约束为，
\[ t_ww_tH_t=T_t \]

消费者最大化问题为，
\begin{align*}
\max_{c_t,h_t}\qquad&\sum_{t=0}^\infty u(c_t^i,h_t^i)	\\
s.t.\qquad& c_t^i=(1-t_w)w_th_t^i+r_tk_t^i+T_t-I_t^i\\
&w_t=f_h(K_t,H_t)\\
&r_t=f_k(K_t,H_t)\\
&k_{t+1}^i=(1-\delta)k_t^i+I_t^i
\end{align*}
或者将约束条件进行一系列替换而简化得到，
\begin{align}
	\max_{c_t,h_t}\qquad&\sum_{t=0}^\infty \beta^tu(c_t^i,h_t^i)	\\\label{inf_eq_cnt}
	s.t.\qquad& c_t^i=(1-t_w)f_h(K_t,H_t)h_t^i+f_k(K_t,H_t)k_t^i+T_t-[k_{t+1}^i-(1-\delta)k_t^i]
\end{align}
这样，就可以写出拉格朗日函数，
\[ L= \sum_{t=0}^\infty \beta^t\{u(c_t^i,h_t^i)-\lambda_t[(1-t_w)f_h(K_t,H_t)h_t^i+f_k(K_t,H_t)k_t^i+T_t-k_{t+1}^i+(1-\delta)k_t^i-c_t^i]\}\]
然后通过选择$c_t^i, h_t^i,k_{t+1}^i $而得到三个最优化条件如下，
\begin{align*}
	\frac{\partial L}{\partial c_t^i}=&u_c(c_t^i,h_t^i)+\lambda_t=0\\
\frac{\partial L}{\partial h_t^i}=&		u_h(c_t^i,h_t^i)-\lambda_t(1-t_w)f_h(K_t,H_t)=0\\
\frac{\partial L}{\partial k_{t+1}^i}=&\beta^t\lambda_t-\beta^{t+1}\lambda_{t+1}[f_k(K_t,H_t)+1-\delta]=0
\end{align*}
消去$ \lambda $则有，
\begin{align*}
	\frac{u_h(c_t^i,h_t^i)}{u_c(c_t^i,h_t^i)}=&(t_w-1)f_h(K_t,H_t)\\
	\frac{u_c(c_t^i,h_t^i)}{u_c(c_{t+1}^i,h_{t+1}^i)}=&\beta [f_k(K_{t+1},H_{t+1})+1-\delta]
\end{align*}
因为个体一样，我们可以把上述一阶条件直接写成总量形式，
\begin{align}\label{inf_eq_h}
	\frac{u_h(C_t,H_t)}{u_c(C_t,H_t)}=&(t_w-1)f_h(K_t,H_t)\\\label{inf_eq_k}
	\frac{u_c(C_t,H_t)}{u_c(C_{t+1},H_{t+1})}=&\beta [f_k(K_{t+1},H_{t+1})+1-\delta]
\end{align}
再来加总约束条件\eqref{inf_eq_cnt}式，有，
\begin{align}\nonumber
&C_t= (1-t_w)f_h(K_t,H_t)H_t+f_k(K_t,H_t)K_t+T_t-[K_{t+1}-(1-\delta)K_t]\\\label{inf_eq_K}
\Longrightarrow 	&K_{t+1}=(1-\delta)K_t+f(K_t,H_t)-C_t\qquad\text{用到了零利润条件的生产函数}
\end{align}
注意\eqref{inf_eq_h},\eqref{inf_eq_k},\eqref{inf_eq_K}这三个一阶条件，与前面没有扭曲时的竞争均衡的三个条件进行比较，会发现其他两个条件是一致的，但\eqref{inf_eq_h}式与前述对应条件是不同的。有了税收以后，个体的效用会降低。

\chapter{确定性递归模型}
尽管变分法可以得到稳态，但有个时候，我不仅禁感兴趣稳态，我们也想知道稳态之前的运动路径，这就需要另外一种递归的方法，动态最优控制理论。
\section{状态和控制}
对一个一个简单的单人RC模型，鲁滨孙需要最大化其贴现后的毕生效用函数，
\begin{align*}
max\; & \sum_{t=0}^\infty \beta^tu(c_{t})\\
s.t.\;& k_{t+1}=(1-\delta )k_t+i_t\\
& y_t=f(k_t)=c_t+i_t,\;\;\;t=1,\cdots,\infty
\end{align*}

控制变量的选择往往对解题难易产生重大影响。但要记住，无论选择谁作为控制变量，都必须有足够的市场或者预算约束使得可以确定时期t其他相关变量的取值。
本例中，$c_t $或者$k
_{t+1} $都可以作为控制变量。
当选择好了$c_t $，根据第二个约束，就选择好了$i_t $，再根据第一个约束就选择好了$k
_{t+1} $ 。
当选择好了$k
_{t+1} $，根据第一个约束就选好了$i_t $	，再根据第二个约束就选好了$c_t $。数学上，这就表现为
\[c_t=f(k_t)+(1-\delta)k_t-k_{t+1}\]

于是目标函数可以写为，
\[max\; \sum_{t=0}^{\infty}\beta^tu[f(k_t)+(1-\delta)k_t-k_{t+1}]\]

\section{值函数}
如果现在我们来控制$k_{t+1} $，那么效用的值函数为，
\begin{align}\nonumber
V(k_t)=\max_{\{k_s\}_{s=t+1}^\infty}&\{u[f(k_t)+(1-\delta)k_t-k_{t+1}]+\beta\cdot u[f(k_{t+1})+(1-\delta)k_{t+1}-k_{t+2}]+\\\nonumber
&\beta^2\cdot u[f(k_{t+2})+(1-\delta)k_{t+2}-k_{t+3}]+\cdots\}\\\label{abc_eq1}
=\max_{\{k_s\}_{s=t+1}^\infty}&\sum_{i=0}^\infty \{\beta^iu[f(k_{t+i})+(1-\delta)k_{t+i}-k_{t+i+1}]\}
\end{align}

用$V(k_t) $表示贴现后的效用值，是为了强调它是初始资本存量的函数。稍微改变下标，将其往前推进一期，有，
\begin{equation}\label{rbc_v1}
V(k_{t+1})=\max_{\{k_s\}_{s=t+2}^\infty}\sum_{i=0}^\infty \beta^iu(f(k_{t+i+1})+(1-\delta)k_{t+1+i}-k_{t+i+2})	
\end{equation}

对于\eqref{abc_eq1}式可以将$t$期的表达剥离出来，写成一个$t$期与$t$以后时期的一个表达如下，
\begin{align}\nonumber
V(k_t)= & \max_{k_{t+1}}[u(f(k_t)+(1-\delta)k_t-k_{t+1})\\\label{rbc_v}
& +\max_{\{k_s\}_{s=t+2}^\infty}\sum_{i=1}^\infty \beta^iu(f(k_{t+i})+(1-\delta)k_{t+i}-k_{t+i+1})]
\end{align}

上式第二部分可以写为，
\begin{align*}
	\max_{\{k_s\}_{s=t+2}^\infty}&\sum_{i=1}^\infty \beta^iu(f(k_{t+i})+(1-\delta)k_{t+i}-k_{t+i+1})=	\max_{\{k_s\}_{s=t+2}^\infty}\{\beta\cdot u[f(k_{t+1})+(1-\delta)k_{t+1}-k_{t+2}]+\\
	&\beta^2\cdot u[f(k_{t+2})+(1-\delta)k_{t+2}-k_{t+3}]+\beta^3\cdot u[f(k_{t+3})+(1-\delta)k_{t+3}-k_{t+4}]+\cdots\}\\
	&=\max_{\{k_s\}_{s=t+2}^\infty}\{\beta\{u[f(k_{t+1})+(1-\delta)k_{t+1}-k_{t+2}]+\beta\cdot u[f(k_{t+2})+(1-\delta)k_{t+2}-k_{t+3}]+\\
	&\beta^2\cdot u[f(k_{t+3})+(1-\delta)k_{t+3}-k_{t+4}]+\cdots\}\}\\
	&=\max_{\{k_s\}_{s=t+2}^\infty}\{\beta\sum_{i=0}^\infty\beta^iu[f(k_{t+1+i})+(1-\delta)k_{t+1+i}-k_{t+i+2}]\}
\end{align*}
因此，结合式\eqref{rbc_v1}，可以把式\eqref{rbc_v}重新写为，
\begin{align*}
V(k_t)= & \max_{k_{t+1}}[u(f(k_t)+(1-\delta)k_t-k_{t+1})\\
& +\beta \max_{\{k_s\}_{s=t+2}^\infty}\sum_{i=0}^\infty \beta^iu(f(k_{t+1+i})+(1-\delta)k_{t+1+i}-k_{t+i+2})]
\end{align*}
即，
\begin{equation}\label{abc_eq2}
V(k_t)=  \max_{k_{t+1}}[u(f(k_t)+(1-\delta)k_t-k_{t+1})+\beta V(k_{t+1})]
\end{equation}

形如\eqref{abc_eq2}式的方程就称为Bellman方程。它表达的问题与\eqref{abc_eq1}式相同，只是它是一个递归形式。它将一个多时期的动态优化变成了通过选择$k_{t+1} $的单期的优化问题。对上式一阶求导并令之为0，可得到，
\begin{equation}\label{abc_eq3}
-u'[f(k_t)-k_{t+1}+(1-\delta )k_t]+\beta V'(k_{t+1}))=0
\end{equation}
实际上，问题并没有变得更简单，因为$V'(k_{t+1}) $未知。然而包络定理表明值函数对某参数的导数等于拉格朗日函数对该参数的偏导在最优处的取值。那么观察\eqref{abc_eq2}式，在那里，$k_{t+1} $是控制变量，$k_t $是一个参数，于是根据包络定理，
\[
V'(k_t)=u'[f(k_t)-k_{t+1}+(1-\delta )k_t]
(f'(k_t)+(1-\delta))\]
即
\[
V'(k_{t+1})=u'(c_{t+1})(f'(k_{t+1})+(1-\delta))
\]
代入\eqref{abc_eq3}式，可以得到
\[
\frac{u'(c_t)}{u'(c_{t+1})}=\beta(f'(k_{t+1})+(1-\delta))
\]
稳态时，$c_t=c_{t+1} $，故稳态的欧拉方程为，
\[
\frac{1}{\beta}-(1-\delta)=f'(\bar k)
\]
\section{一般形式}
令$x_t $为$t$时期状态变量所组成的向量，$ y_t $为$t$时期控制变量所组成的向量。$F(x_t,y_t) $为$t$时期需要最大化的目标函数值。在上例中，给定初始状态变量$ x_t $的值，要求解的问题是值函数，
\[
V(x_t)=\max_{\{y_s\}_{s=t}^\infty}\sum_{s=t}^\infty \beta^{s-t}F(x_s,y_s)
\]
其预算约束集为所有$s\ge t $时期的
\[
x_{s+1}=G(x_s,y_s)
\]

在满足$s\ge t $各个时期，目标函数$F$和约束$G$形式都是相同的。那么根据前述讲解，值函数就可以写成递归形式，即Bellman方程，
\[
V(x_t)=\max_{y_t}[F(x_s,y_s)+\beta V(x_{t+1})]
\]

满足约束，
\[
x_{s+1}=G(x_s,y_s)
\]

或者将约束代入目标函数，可得
\begin{equation}\label{abc_eq4}
V(x_t)=\max_{y_t}[F(x_t,y_t)+\beta V(G(x_t,y_t))]
\end{equation}

求解就可以得到控制变量关于t时期状态变量的函数
\[
y_t=H(x_t)
\]

此即政策函数。如果假设我们已经得到了政策函数，那么将政策函数代入\eqref{abc_eq4}式，就有
\[
V(x_t)=F(x_t,y_t)+\beta V[G(x_t,H(x_t))]
\]

上式不再有最大化符号，因为这种优化已经隐含在政策函数中了。
所以，关键在于求得政策函数。对\eqref{abc_eq4}式按照控制变量得到一阶条件，
\begin{equation}\label{abc_eq5}
F_y(x_t,y_t)+\beta V'(G(x_t,y_t))G_y(x_t,y_t)=0
\end{equation}

其中，$F_y $是目标函数对控制变量的偏导，$V'(G(x_t,y_t)) $是值函数对$t+1$期状态变量的偏导，$G_y $是预算约束对控制变量的偏导。关键在于$V'(G(x_t,y_t)) $，一般是用包络定理来求。包络定理表明，
\[
V'(x_t)=F_x(x_t,y_t)+\beta V'(G(x_t,y_t))G_x(x_t,y_t)
\]

针对这样的表达式，可以通过对控制变量的艺术选择，使得$G_x(x_t,y_t)=0 $，则上式可以简化为，
\[
V'(x_t)=F_x(x_t,y_t)
\]

一阶条件\eqref{abc_eq5}式就可以简化为，
\begin{equation}\label{abc_eq6}
F_y(x_t,y_t)+\beta F_x[G(x_t,y_t),y_{t+1}]G_y(x_t,y_t)=0
\end{equation}

现在就可以分情况来讨论，
\begin{enumerate}
	\item 如果$F_x[G(x_t,y_t),y_{t+1}] $与$y_{t+1} $无关，那么\eqref{abc_eq6}式就给出了政策函数的隐函数形式，于是就可以得到值函数的隐函数形式。
	\item 如果$F_x[G(x_t,y_t),y_{t+1}] $与$y_{t+1} $相关，那么可以按照$y_{t+1}=y_t $来求出稳态情况。
\end{enumerate}

如果非常遗憾地，对控制变量的选择无法达到$G_x(x_t,y_t)=0 $，则需要数值解来得到值函数。数值解的步骤如下，
\begin{enumerate}
	\item 给出一个值函数的初始猜测$V_0(x_t) $，初值是什么不重要，所以有时干脆令之为0。
\item	使用$V_1(x_t)=\max_{y_t}[F(x_t,y_t)+\beta V_0(G(x_t,y_t))] $来得到新的值函数$V_1(x_t) $。在$x_t$的定义域内给出一系列的点，然后针对每个点运用计算机最大化算法得到数值。以此类推，在某些条件满足的前提下，该序列会收敛到值函数$V(x_t) $。	
\end{enumerate}

实际上，该过程中，政策函数也被描绘。每一个给定的$ x_t $所对应的序列$y_t$的极限点，精确的刻画了政策函数$y_t=H(x_t) $。	

\section{回到实例经济}
回到本章开始的最优化问题，其目标函数为
\[
F(x_t,y_t)=u(f(k_t)-k_{t+1}+(1-\delta)k_t)
\]

所谓预算约束的目的就是要通过$t$期的控制和状态变量来确定$t+1$期的状态变量。本例中，十分特殊，$t$期的控制变量就是$t+1$的状态变量，因此，不需要具体的函数形式$G$，只要通过$t$期状态和控制变量，就已经知道了$t+1$期的状态变量。

如果选择$ c_t $作为控制变量，那么需要重新定义目标函数和预算约束。此时目标函数为，
\[
F(x_t,y_t)=u(c_t)
\]

预算约束为
\[
k_{t+1}=f(k_t)+(1-\delta)k_t-c_t
\]

对应的Bellman方程为，
\[
V(k_t)=max_{c_t}\; [u(c_t)+\beta V(f(k_t)+(1-\delta)k_t-c_t)]
\]

但此时使用包络定理，得到
\[
\frac{\partial G(x_t,y_t)}{\partial x_t}=f'(k_t)+(1-\delta)
\]

于是，
\begin{align*}
V'(x_t) & =F_x(x_t,y_t)+\beta V'(G(x_t,y_t))G_x(x_t,y_t)\\
& = \beta V'(G(x_t,y_t))(f'(k_t)+(1-\delta))\\
\end{align*}

值函数又被转换成了值函数的导数及其他项的表达式，没有任何改进。为了避免这种情况，最好就是尽可能多的把模型都放进目标函数，\textbf{让预算约束中不再含有$t$时期的状态变量}。
\section{数值求解}
\subsection{值函数迭代}
考虑具体的函数形式，生产函数为，
\[
f(k_t)=k_t^{\theta}
\]

其中$0<\theta<1 $，效用函数为
\[
u(c_t)=ln(c_t)
\]

对应的Bellman方程为，
\[
V(k_t)=max_{k_{t+1}}\; [ln(k_t^\theta-k_{t+1}+(1-\delta)k_t)+\beta V(k_{t+1})]
\]

为使可以使用递归方法逼近$V(\cdot) $，需要各参数值，然后进行迭代。看程序\lstinline|val_it.R|（(它有子函数\lstinline|valfun.R|)会比较清晰。
\subsection{欧拉方程迭代}
用欧拉方程迭代，会比值函数迭代速度要快好几倍。该最优规划最终的一阶方程组如下，第一个方程是资源约束方程，第二个方程是欧拉方程。

\begin{align}\label{abc_eq_sr}
	\underbrace{k_{t+1}}_{a_t}=&\underbrace{f(k_t)+(1-\delta)k_t}_{m_t}-c_t	\\\label{abc_eq_eu}
	u'(c_t) =&\beta u'(c_{t+1})\cdot [f'(k_{t+1})+(1-\delta)]
\end{align}

具体算法如下(代码见\lstinline|EGM.R|)：
\begin{enumerate}
	\item 格点化$ k_t $，猜测一个与$ k_t $格点数目一样的$ c_t $序列，形成数值对序列$ \{k_t,c_t^0\} $
	\item 基于\eqref{abc_eq_sr}式，利用$k_t, c_t^0 $,计算得到$ k_{t+1}^0 $
	\item 利用数值对$ \{k_t,c_t^0\} $和$ k_{t+1}^0 $插值得到$ c_{t+1}^0 $
	\item 再利用\eqref{abc_eq_eu}式反算出$ c_t^1 $，更新$ c_t^0 $
	\item 重复第2步至第4步，直至收敛
\end{enumerate}

我把这个思路总结一下：
\begin{note}
	首先把状态变量$ k_t $格点化，为每个格点随意猜测一个控制变量序列$ c_t $；然后，基于这个状态变量的格点以及控制变量的猜测序列，利用差分方程组，得到下一期的状态变量$ k_{t+1} $；然后，通过插值得到新的$ c_{t+1} $，再次利用差分方程组又可以反算出上一期的控制变量，即$ c_t $。也就是说，我们从$ c_t $的猜测出发，利用差分方程组又算出来一个$ c_t $，我们需要这个$ c_t $是一样的，那么只要往复迭代，就能最终找到真正的$ c_t $。我们把这个过程简要总结在图\ref{abc_fig_alg}中。
	\end{note}

\begin{figure}[H]\centering
\begin{tikzpicture}[inner sep=1mm,
	place/.style ={rectangle,fill =black!20, draw = black, minimum size = 1cm},
	arrow1/.style = {draw = blue!50,-{Latex[length = 2mm, width = 1mm]}},]
	\node[place] (A) {$ (k_t,c_t^0) $};
	\node[place] (B) [right = 3cm of A]{$ k_{t+1}^0 $};
	\node[place] (C) [right = 2cm of B]{$ c_{t+1}^0 $};
	\node[place] (D) [right = 2cm of C]{$ c_{t}^1 $};
	\draw[arrow1] (A) to node[above]{资源约束方程} (B);
	\draw[arrow1] (B) to node[above]{插值} (C);
	\draw[arrow1] (C) to node[above]{欧拉方程} (D);
	\draw[arrow1] (D) to[out=90,in =90] node[above]{迭代直至收敛} (A);
\end{tikzpicture}
\caption{算法图示}\label{abc_fig_alg}
\end{figure}

几点评论：
\begin{itemize}
	\item 我们最终获得数值对$ \{k_t,c_t\} $，实际上就是政策函数
\end{itemize}

\subsection{内生格点法}
代码见\lstinline|EGM.R|。
\begin{enumerate}
	\item 格点化$ a_t $，猜测一个数值对序列$ \{m_t^0,c_t^0\} $
	\item 因为$ a_t $就是$ k_{t+1} $，因此，利用$ a_t $可以得到$ m_{t+1}^{0} $，然后基于数值对$ \{m_t^0,c_t^0\} $插值得到$ c_{t+1}^0 $
	\item 利用$ a_t,c_{t+1}^0 $和\eqref{abc_eq_eu}式得到$ c_t^1 $，从而更新了$ c_t^0 $
	\item 利用$ m_t^1=a_t+c_t^1 $更新$ m_t^0 $
	\item 重复第2步至第4步，直至收敛
\end{enumerate}

\textbf{可以看到标准的解法和内生格点法真正的区别在于因为$ a_t $的设置，使得期望的计算大为简化，要知道如果\eqref{abc_eq_eu}式中涉及到t+1的脚标都是要计算期望的\footnote{我这个例子用的是确定性模型，随机模型是必须要计算期望的。}，标准解法必须计算所有格点上的期望，而内生格点法只需要在插值$ c_{t+1} $的时候计算一下。}

\section{相图与政策函数}
相图有助于我们理解差分方程的运行规律。所谓相图只适用于$ 2\times 2 $的微分方程组，它是以微分方程组的内生变量为横轴和纵轴绘制的定性图形。

针对前面的例子，具体操作如下，将\eqref{abc_eq_sr}式和\eqref{abc_eq_eu}式重新写成
\begin{align*}
	k_{t+1}=&k_t^\theta+(1-\delta)k_t-c_t	\\
	\frac{1}{c_t} =&\beta \frac{1}{c_{t+1}}\cdot [\theta k_{t+1}^{\theta-1}+(1-\delta)]
\end{align*}
然后做一个移项，
\begin{align}\label{abc_eq_sr2}
	k_{t+1}-k_t=&k_t^\theta-\delta k_t-c_t	\\\label{abc_eq_eu2}
	c_{t+1}-c_t =&\beta c_t\cdot (\theta k_{t+1}^{\theta-1}+1-\delta)-c_t
\end{align}
然后去掉方程组的时间脚标，这样方程组左边的差分项为0，这样就可以得到两个非线性方程，它们对应的是单变量的稳态。每个非线性方程，我们可以在坐标中绘制一根曲线。
\begin{figure}[H]
\includegraphics[scale=0.7]{DynamicMacroeconomicsRCode/data-raw/phase.pdf}
\caption{相图(见\lstinline|EGM.R|)}
\end{figure}

\begin{itemize}
	\item 曲线就是第一个方程，垂线就是第二个方程
	\item 带箭头的就是政策函数，也是我们通常说的鞍形路径
	\item 通过这个图，我们可以知道图像中任何一点会如何运动
	\begin{itemize}
		\item 比如对于蓝色曲线上下，变量如何运动？因为曲线上的点意味着$ k $是不变的，那么假如曲线上某一点稍微偏离该曲线比如往下一点点，从\eqref{abc_eq_sr2}式可以知道此时$ k $会变大，这就表明曲线下面就意味着$ k $向右边运动。
		\item 类似地，现在垂线上的点意味着$ c $是不变的，那么假如垂线上一点往右边变化一点点，从\eqref{abc_eq_eu2}式可以知道$ c_t $会变小，这就表明垂线右边就意味着$ c $往下运动，左边$ c $往上运动。
	\end{itemize}
\end{itemize}
\section{一个可变劳动的例子}
个人最优规划可以写为，
\begin{align*}
	\max\sum_{i=0}^\infty \qquad&\beta^iu(c_{t+1},h_{t+1})\\
	s.t.\qquad& k_{t+1}=(1-\delta)k_t+i_t\\
	& f(k_t,h_t)=c_t+i_t\\
	& h_t\le 1
\end{align*}

写成贝尔曼方程形式，
\[ V(k_t)=\max_{k_{t+1},h_t}\{u[f(k_t,h_t)+(1-\delta)k_t-k_{t+1},h_t]+\beta V(k_{t+1})\} \]

实际上，当你一旦有了贝尔曼方程就可以数值求解了。一个具体使用值函数迭代的例子见代码\lstinline|valfun_varlab.R|。

接下来，看看它的一阶条件。注意到它的约束$ k_{t+1}=G(k_{t+1})=k_{t+1}$，因此$ G_x(x_t,y_t)=0 $。此时两个控制变量$ k_{t+1},h_t $，因此有两个一阶条件，即
\begin{align*}
	\frac{\partial V(k_t)}{\partial h_t}=&u_c(f(k_t,h_t)+(1-\delta)k_t-k_{t+1},h_t)\cdot f_h(k_t,h_t)+u_h(f(k_t,h_t)+(1-\delta)k_t-k_{t+1},h_t)=0\\
	\frac{\partial V(k_t)}{\partial k_{t+1}}=&-u_c(f(k_t,h_t)+(1-\delta)k_t-k_{t+1},h_t)+\beta V'(k_{t+1})
\end{align*}
而包络定理意味着，
\[ V'(k_t)=u_c[f(k_t,h_t)+(1-\delta)k_t-k_{t+1},h_t]\cdot [f_{k}(k_t,h_t)+1-\delta] \]
利用上面三个等式，同样可以得到，
\begin{align}\label{rec_eq_k}
	\frac{u_h(c_t,h_t)}{u_c(c_t,h_t)}&=-f_h(k_t,h_t)\\\label{rec_eq_c}
	\frac{u_c(c_t,h_t)}{u_c(c_{t+1},h_{t+1})}&=\beta[f_k(k_{t+1},h_{t+1})+(1-\delta)]
\end{align}
这两个式子和前面第\ref{inf}章中的\eqref{vl7}式和\eqref{vl8}式是一样的。如果此时设置效用函数和生产函数为，
\begin{align*}
	u(c_t,h_t)=&\ln c_t+A\ln(1-h_t)\\
	f(k_t,h_t)=&k_t^\theta h_t^{1-\theta}
\end{align*}
那么模型的一阶条件和约束条件可以写为，
\begin{align}\label{rec_eq_h}
	\frac{Ac_t}{1-h_t}=&(1-\theta)k_t^\theta h_t^{-\theta}\\\label{rec_eq_cc}
	\frac{c_{t+1}}{c_t}=&\beta (\theta k_t^{\theta-1}h_t^{1-\theta}+1-\delta)\\\label{rec_eq_knext}
	k_{t+1}=&(1-\delta)k_t+k_t^\theta h_t^{1-\theta}-c_t
\end{align}

对于这个非线性差分方程组也可以使用内生格点法(见代码\lstinline|val_ite_varlab.R|)：
\begin{enumerate}
	\item 格点化$ k_t $，并猜测控制变量$ c_t,h_t $，形成数值对$ (k_t,c_t^0,h_t^0) $
	\item 利用数值对$ (k_t,c_t^0,h_t^0) $以及\eqref{rec_eq_knext}式得到$ k_{t+1}^0 $
	\item 利用数值对$ (k_t,c_t^0) $插值获得$ c_{t+1}^0 $
	\item 利用\eqref{rec_eq_cc}式和$ c_{t+1}^0,k_t,h_t^0 $获得更新后的$ c_t^1 $
	\item 利用\eqref{rec_eq_h}式和$ c_t^1,k_t $获得更新后的$ h_t^1 $
	\item 重复上述步骤直至收敛
\end{enumerate}

几个值得注意的事项：
\begin{itemize}
	\item 注意到\eqref{rec_eq_h}式涉及到非线性方程求根，在R语言中使用\lstinline|rootSolve::unitroot.all|函数，不要用\lstinline|unitroot|。因为前者会在你给定的区间主动再次细分区域，找到全部的根，算法更容易成功。
	\item 在计算的过程中，对$ c_t,h_t $的猜测实际上也不是任意的，猜测的不好会导致方程无法收敛。比如在该方程中，如果$ c_t $猜测的过大，会使得利用\eqref{rec_eq_knext}式计算$ k_{t+1} $时变成负数\footnote{该类方程报错一般都在此处，要特别注意此处。实际上也可以通过先算绝对值再加负号的方法迂回得到。比如\lstinline|(-0.3)^(1/3)|是算不出来的，但是可以\lstinline|-(0.3)^(1/3)|。}，负数在此处也不是一个有意义的值，不通过处理，方程会发散、报错。
\end{itemize}


\chapter{随机性递归模型}\label{sto}
求导符号可以穿过期望算子。

\section{一个简单的随机增长模型}
相对于以前的RC模型，我们仅修改生产函数为，
\[ y_t=A^tf(k_t) \]

其中$ A^t $是技术水平，但是它随机取两个值，以概率$ p_1 $取$ A_1 $，以概率$ p_2 $取$ A_2 $。可以把它看成农业生产中的天气变化。然后资本存量增长方程也有细微变化，
\[ k_{t+1}=A^tf(k_t)+(1-\delta)k_t-c_t \]
个体最大化贴现后的期望效用$ E_0\sum_{t=0}^\infty \beta^t u(c_t) $，那么在初期0期，值函数及其所受约束可以写成，
\begin{align*}
	V(k_0,A^0)=&\max_{\{c_t\}_{t=0}^\infty}E_0\sum_{t=0}^\infty\beta^tu(c_t)\\
	s.t.\qquad k_{t+1}=&A^tf(k_t)+(1-\delta)k_t-c_t
\end{align*}

要注意，现在的值函数与以前的不同在于，多了一个状态变量$ A^0 $，这个东西会给我们数值求解带来一定的麻烦。现在按照前面的过程，来把它写成递归的形式为，
\begin{align*}
	V(k_0,A^0)=&\max_{c_0} [u(c_0)+\beta E_0V(k_1,A^1)]\\
	s.t.\qquad k_{1}=&A^0f(k_0)+(1-\delta)k_0-c_0
\end{align*}
或者，我们可以写成更加通用的形式，不是第0期，而是第$ t $期，那么，
\begin{align*}
	V(k_t,A^t)=&\max_{c_t} [u(c_t)+\beta E_tV(k_{t+1},A^{t+1})]\\
	s.t.\qquad k_{t+1}=&A^{t}f(k_t)+(1-\delta)k_t-c_t
\end{align*}

如果选择$ k_{t+1} $为控制变量，那么我们要求解的问题可以写成，
\begin{align}\label{sto_eq_v}
	V(k_t,A^t)=&\max_{k_{t+1}} [u(A^{t}f(k_t)+(1-\delta)k_t-k_{t+1})+\beta E_tV(k_{t+1},A^{t+1})]\\
	s.t.\qquad k_{t+1}=&G(x_t,y_t)=k_{t+1}
\end{align}

此时，数值求解值函数最大的问题是维数问题。回忆前面我们数值求解时，为了画出值函数，我们在状态变量$ k_t $的连续定义域上选取了从以0.06为格子选取了一系列的点，然后基于这些点，对目标函数进行了最优化而得到值函数，在这个过程中，我们还用线性插值来计算这些离散点之间的目标函数值。当时我们就仅仅计算每个格点处的值就可以，我们的代码当时是100个离散点，或者说维度$ dim(x_t)=100 $。

现在有了随机状态变量，如果该随机状态变量只取2个值，那么计算值函数时，每个状态变量$ k_t $的离散点上，必须要计算其对应随机状态的两个值，也就是说此时的计算维度是$ dim(x_t)\cdot dim(A^t) $。

我们来试一下，针对前面的例子如何得到值函数。
\section{简单经济的值函数}
在数值模拟中，要将\eqref{sto_eq_v}式中的期望符号换成具体可计算的形式。也就是说针对$ A^t $以概率$ p_1,p_2 $的两个取值$ A_1,A_2 $，\eqref{sto_eq_v}式实际上要写成两个式子，
\begin{align}\label{sto_eq_A1}
	V(k_t,A_1)=&\max_{k_{t+1}} [u(A_1f(k_t)+(1-\delta)k_t-k_{t+1})+\beta [p_1V(k_{t+1},A_1)+p_2V(k_{t+1},A_2)]]\\\label{sto_eq_A2}
		V(k_t,A_2)=&\max_{k_{t+1}} [u(A_2f(k_t)+(1-\delta)k_t-k_{t+1})+\beta [p_1V(k_{t+1},A_1)+p_2V(k_{t+1},A_2)]]
\end{align}

回忆值函数迭代过程，首先是要给定初值$ V_0(k_t,A_1),V_0(k_t,A_2) $，然后依据\eqref{sto_eq_A1}和\eqref{sto_eq_A2}式，就可以迭代得到下一轮的$ V_1(k_t,A_1),V_1(k_t,A_2) $，以此类推，直到收敛。
\begin{figure}[H]
\includegraphics[scale=0.7]{val_ite_sto.png}
\includegraphics[scale=0.7]{policy_sto.png}
\caption{值函数迭代和政策函数}
来自源码\lstinline|val_ite_sto.R|和\lstinline|valfun_sto.R|。
\end{figure}

一旦有了政策函数，实际上就可以模拟经济中比如资本存量随时间的演化过程。比如给定一个初始值$ k_0=1 $，因为政策函数是把当期状态变量和当期控制变量联系起来的函数。在我们这个例子里，当期状态变量是$ k_t $，当期的控制变量正好是$ k_{t+1} $，这样根据政策函数，给定$ k_0 $，就能得到$ k_1 $；给定$ k_1 $，就能得到$ k_2 $等等。当然，因为我们的政策函数实际上是一系列的点，毫无疑问，我们又需要插值。我们在代码\lstinline|val_ite_sto.R|的最后部分绘制了该经济体$ k_t $从1开始的演变路径。注意，因为不同状态下，有两个政策函数，我们设置了一个均匀分布的随机数发生器，随机数小于等于0.8时，使用状态1的政策函数，否则，使用状态2的政策函数。
\begin{figure}[H]
\includegraphics[scale=0.7]{simu_k.png}
\caption{模拟出来的时间路径}
\end{figure}

对此模型，也可以使用差分方程组迭代，而不用值函数迭代。写出一阶非线性差分方程组，
\begin{align*}
	u'(c_t)=&\beta E_t[u'(c_{t+1})(f'(k_{t+1})+1-\delta)]\\
	k_{t+1}=&A^{t}f(k_t)+(1-\delta)k_t-c_t
\end{align*}
对于具体的函数形式，$ f(k_t)=A^tk_t^\theta,u=\ln c_t $，上述式子可以写成，
\begin{align}\label{sto_eq_cc}
	\frac{1}{c_t}=&\beta E_t\left[\frac{1}{c_{t+1}}(A^{t+1}\theta k^{\theta-1}_{t+1}+1-\delta)\right]\\\label{sto_eq_ktp1}
	k_{t+1}=&A^{t}k_t^\theta+(1-\delta)k_t-c_t
\end{align}

使用差分方程迭代时，要注意期望的计算，这里涉及到2期的概率组合。具体而言，在格点化$ k_t $后，利用\eqref{sto_eq_ktp1}式计算$ k_{t+1} $时，涉及到$ A^t $的出现概率，其可为1.75或者0.75，然后在利用\eqref{sto_eq_cc}式计算期望时，涉及到$ A^{t+1} $，其可为1.75或0.75，因此，为得到\eqref{sto_eq_cc}式的期望，我们涉及到$ 2\times 2=4 $个组合，即$ \{A_1A_1,A_1A_2,A_2A_1,A_2A_2\} $。这部分关键代码我们书写如下所示(完整代码见\lstinline|val_ite_sto.R|)，因为涉及到两个时期的概率，所以我们从第5行代码开始用了2个循环来计算这个期望。
\begin{lstlisting}[language=R,numbers=left,firstnumber=1,mathescape=false]
prob <- data.frame(t1 = c(A1,A1,A2,A2),t2 = c(A1,A2,A1,A2), p = c(p1*p1,p1*p2,p2*p1,p2*p2))
for (i in 1:50) {
	# 计算期望
	EU <- numeric(length(k0))
	for (At1 in c(A1,A2)){# 第t期概率
		knext <- At1*k0^theta + (1-delta)*k0-cc[,ncol(cc)]
		cnext <- interp1(k0,cc[,ncol(cc)],knext,extrap = T)
		for (At2 in c(A1,A2)) {# 第t+1期概率
			EU <- EU + ((At2*theta*knext^(theta-1)+1-delta)/cnext) * prob$p[prob$t1 == At1 & prob$t2 == At2]
		}
	}
	cc <- cbind(cc,1/(beta*EU))
}
\end{lstlisting}


\section{马尔可夫链}
现在我们将随机过程更加复杂一点。

Markov链要记住，它包含三个要素：一是自然状态可能值的集合，比如$ A^t=\{1.75,0.75\} $。二是转移矩阵$ P $,它的元素$ p_{ij} $表达的是上期在状态$ i $下，本期状态取$ j $的概率，很显然，这是一个条件概率。比如，
\[ P=\begin{bmatrix}
	p_{11}&p_{12}\\
		p_{21}&p_{22}
\end{bmatrix}=\begin{bmatrix}
0.9&0.1\\0.4&0.6
\end{bmatrix} \]

三是时期$ t-1 $的初始自然状态，它决定了我们要从$ P $的哪一行开始。

尽管$ P $给出了条件概率矩阵，但我们同样感兴趣无条件概率(边缘概率)。如果经济运行时间够长，无条件概率实际上与初始状态无关，且存在唯一不变的无条件概率。我们不作数学证明，而是从数值上来验算这个收敛如何发生。

回忆计量课中曾经讲过的边缘概率的算法，我们现在已知第1期的概率分布$ p_0=(0.36,0.64)' $来算第2期的边缘概率分布。首先，我们需要知道联合概率，而联合概率有公式$ P(AB)=P(A|B)P(B) $，因此可以得到联合概率分布表\ref{sto_tb1}。然后，就可以计算边缘概率。你很快会发现，第2期的边缘分布就是$ p_0P $，第3期的边缘分布就是$ p_0P ^2$，那么当时期趋于无穷大时，$ p_0P^n\rightarrow P^\infty $，边缘分布与初始状态无关。
\begin{table}[H]\centering
	\caption{联合概率分布表}\label{sto_tb1}
\begin{tabular}{lcc}\hline
& 第2期处于状态1&第2期处于状态2\\\hline
第1期处于状态1&$ p_{11}\cdot0.36 $&$ p_{12}\cdot 0.36 $\\
第1期处于状态2&$ p_{21}\cdot0.64 $&$ p_{22}\cdot0.64 $\\
边缘分布&$p_{11}\cdot0.36+p_{21}\cdot0.64 $&$ p_{12}\cdot 0.36+p_{22}\cdot0.64  $\\\hline
\end{tabular}
\end{table}


对于马尔可夫过程的模型，值函数略有不同，
\[ V(k_t,A^t)=\max_{k_{t+1}} [u(A^{t}f(k_t)+(1-\delta)k_t-k_{t+1})+\beta E_t[V(k_{t+1},A^{t+1})|A^t]] \]

主要区别在于此时的期望是条件期望$E_t[V(k_{t+1},A^{t+1})|A^t]$，因此，具体计算时就要使用条件期望$ p_{ij} $而不是先前的$ p_1,p_2 $。
\begin{align}
	V(k_t,A_1)=&\max_{k_{t+1}} [u(A_1f(k_t)+(1-\delta)k_t-k_{t+1})+\beta [p_{11}V(k_{t+1},A_1)+p_{12}V(k_{t+1},A_2)]]\\
	V(k_t,A_2)=&\max_{k_{t+1}} [u(A_2f(k_t)+(1-\delta)k_t-k_{t+1})+\beta [p_{21}V(k_{t+1},A_1)+p_{22}V(k_{t+1},A_2)]]
\end{align}

依据上述表达，我们可以重新计算值函数和政策函数，还可以继续模拟资本存量的时间路径如下：
\begin{figure}[H]
\includegraphics[scale=0.7]{simu_k_mkv.png}
\caption{马尔可夫过程的资本存量时间路径演化}
\end{figure}

注意模拟马尔可夫过程时，要设定初始状态。具体参见源代码\lstinline|val_ite_mkv.R|。

\chapter{Hansen 的RBC模型}\label{sec1}
第\ref{sto}章讲述的随机化，在维数较小的时候是可行的，一旦维数变大，计算将十分巨大以致根本不可行，需要一些其他的技术做出处理。
\section{Hansen的基本模型}
一个RC经济体，其贴现后效用函数如下，
\[max\; \sum_{t=0}^\infty \beta^tu(c_t,l_t)\]

其中，$c_t $是时期$t$的消费量，$l_t $时期$t$的闲暇且$l_t=1-h_t $，$h_t $是时期$t$的劳动。效用函数的具体形式如下，
\[u(c_t,1-h_t)=lnc_t+Aln(1-h_t)\]

其中$A>0 $。生产函数为包含随机技术的柯布道格拉斯形式，
\[f(\lambda_t,k_t,h_t)=\lambda_t k_t^\theta h_t^{1-\theta}\]

其中$\lambda_t $服从如下随机过程，
\[
\lambda_{t+1}=\gamma \lambda_t +\varepsilon_{t+1}
\]

其中，$0<\gamma<1 $，$\varepsilon_t $独立同分布，正值，上有界，均值为$1-\gamma $。\textbf{这些假设保证了$\lambda_t $的均值为1，且产出不为负值}。资本积累满足下式，
\[k_{t+1}=(1-\delta)k_t+i_t\]

任意时期$t$，有预算约束，
\[f(\lambda_t,k_t,h_t)\ge c_t+i_t\]

\paragraph{上述模型的Bellman方程可以表达为,}
\begin{align*}
V(k_t,\lambda_t)=&\max_{c_t,h_t}\; lnc_t+Aln(1-h_t)+\beta E_t[V(k_{t+1},\lambda_{t+1})|\lambda_t]\\
s.t.\quad & \lambda_t k_t^\theta h_t^{1-\theta}\ge c_t+i_t\\
& \lambda_{t+1}=\gamma \lambda_t+\varepsilon_{t+1}\\
& k_{t+1}=(1-\delta)k_t+i_t
\end{align*}

对于上述最大化表述，我们希望将控制变量$ c_t,h_t $转换成$k_{t+1},h_t $，通过预算约束的代换，可以得到，
\begin{align*}
V(k_t,\lambda_t)=\max_{k_{t+1},h_t}\; & ln[\lambda_t k_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}]+Aln(1-h_t)+\\
& \beta E_t[V(k_{t+1},\lambda_{t+1})|\lambda_t]
\end{align*}

该Bellman方程的一阶条件为，
\[\frac{\partial V(k_t,\lambda_t)}{\partial k_{t+1}}=0=-\frac{1}{\lambda_tk_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}}+\beta E_t[V_k(k_{t+1},\lambda_{t+1})|\lambda_t]\]	
\[\frac{\partial V(k_t,\lambda_t)}{\partial h_{t}}=0=(1-\theta)\frac{1}{\lambda_tk_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}}(\lambda_tk_t^\theta h_t^{-\theta})-A\frac{1}{1-h_t}\]

而包络定理告诉我们，
\[
\frac{\partial V(k_t,\lambda_t)}{\partial k_{t}}=\frac{1}{\lambda_tk_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}}[\theta \lambda_tk_t^{\theta-1} h_t^{1-\theta}+(1-\delta)]
\]
将包络定理代入，一阶条件可以重新书写为，
\begin{equation}\label{eq7}
\frac{1}{\lambda_tk_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}}=\beta E\left[ \frac{\theta \lambda_{t+1}k_{t+1}^{\theta-1} h_{t+1}^{1-\theta}+(1-\delta)}{\lambda_{t+1}k_{t+1}^\theta h_{t+1}^{1-\theta}+(1-\delta)k_{t+1}-k_{t+2}}\left| \lambda_t\right. \right]
\end{equation}

\begin{equation}\label{eq8}
(1-\theta)(1-h_t)(\lambda_t k_t^\theta h_t^{-\theta})=A[\lambda_t k_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}]
\end{equation}
对于可行性约束，有，
\[c_t=\lambda_t k_t^\theta h_t^{1-\theta}+(1-\delta)k_t-k_{t+1}\]
而对于竞争性市场，要素报酬等于其要素边际产品，即，
\[r_t=\theta \lambda_t k_t^{\theta-1} h_t^{1-\theta},\quad w_t=(1-\theta)\lambda_tk_t^\theta h_t^{-\theta}\]
这样，\eqref{eq7}式和\eqref{eq8}式就可以得到一个较为简单的形式，
\[\frac{1}{c_t}=\beta E_t\left [ \frac{r_{t+1}+(1-\delta)}{c_{t+1}}|\lambda_t\right ]\]
\[(1-h_t)w_t=Ac_t\]
令$\bar k=k_t=k_{t+1}=k_{t+2} $，$\bar h=h_t=h_{t+1} $，且稳态时，期望值等于实现值，这样就可以得到稳态的$\bar k $和$\bar h $，
\begin{equation}\label{eq9}
\bar h=\frac{1}{1+\frac{A}{1-\theta}\left[1-\frac{\beta\delta\theta}{1-\beta(1-\delta)} \right]}
\end{equation}
\begin{equation}\label{eq10}
\bar k=\bar h\left[ \frac{\theta \bar \lambda}{\frac{1}{\beta}-(1-\delta)}\right]
\end{equation}
如果把模型写成总量形式（依据福利经济学第二定理），一阶条件\eqref{eq7}和\eqref{eq8}式可以写为，
\[1=\beta E_t\left [ \frac{C_{t}}{C_{t+1}}(r_{t+1}+(1-\delta))\right ]\]
\[(1-H_t)(1-\theta)\frac{Y_t}{H_t}=AC_t\]

预算约束为(不包含存量控制那个约束了)，
\[C_t=Y_t+(1-\delta )K_t-K_{t+1}\]

产品市场上两个约束为，
\[Y_t=\lambda_tK_t^\theta H_t^{1-\theta}\]
\[r_t=\theta \frac{Y_t}{K_t}\]

实际上，我们十分关心稳态之外的方程运动方式。本来第\ref{sto}章的方法也可以用，无非就是不断迭代，算出值函数及政策函数。但问题在于，如果随机变量的维数迅速增大，这种迭代值函数的办法会变得无法负担，我们需要一些线性化的方法。一种就是下一小节要阐述的对数线性化，另一种就是下章要阐述的二次规划。

\section{ 对数线性化方法}
\subsection{对数线性化基础}
该方法分两步，第一步是对函数取对数，第二步就是对对数化后的式子一阶泰勒展开。
本小节内容简单，且今后使用较少。重点看看Uhig的对数线性化方法。
\subsection{Uhlig的对数线性化方法}
如果定义
\[\widetilde X_t=lnX_t-ln\overline X\]

就有
\[X_t=\overline Xe^{\widetilde X_t} \]

于是(一阶泰勒展开得到，p61)，
\begin{align*}
e^{\widetilde X_t+a\widetilde Y_t}&\approx 1+\widetilde X_t+a\widetilde Y_t\\
\widetilde X_t\widetilde Y_t&\approx 0\hspace{2em}\text{高阶的消去而已}\\
E_t[ae^{\widetilde X_{t+1}}]&\approx a+aE_t(\widetilde X_{t+1})	\hspace{2em}\text{第一个法则在期望里面展开而已}\\
\end{align*}

第一个法则是经常使用的。

\section{Hansen模型的对数线性形式}
之所以要对数线性化，是得到稳态之后，想知道稳态之外的运动形式。稳态是提前已经求出了的。那么对于Hansen模型，其一阶条件以及相关预算约束共5个方程可以重新书写如下，
\begin{align*}
&1=\beta E_t\left [ \frac{C_{t}}{C_{t+1}}(r_{t+1}+(1-\delta))\right ]\\
&(1-H_t)(1-\theta)\frac{Y_t}{H_t}= AC_t\\
&C_t=Y_t+(1-\delta )K_t- K_{t+1}\\
&Y_t=\lambda_tK_t^\theta H_t^{1-\theta}\\
&r_t=\theta \frac{Y_t}{K_t}
\end{align*}

现在使用对数线性化方法，首先进行$X_t=\overline Xe^{\widetilde X_t} $ 变换，然后进行$e^{\widetilde X_t+a\widetilde Y_t}\approx 1+\widetilde X_t+a\widetilde Y_t $变换。那么上述5个方程经过变化后，分别为，
\begin{align*}
&0\approx \widetilde C_t-E_t\widetilde C_{t+1}+\beta \overline r E_t\widetilde r_{t+1}\\
&0\approx \widetilde Y_t-\frac{\widetilde H_t}{1-\overline H}-\widetilde C_t\\
&0\approx \overline Y\widetilde Y_t-\overline C\widetilde C_t+\overline K[(1-\delta)\widetilde K_t-\widetilde K_{t+1}]\\
&0\approx \widetilde \lambda_t+\theta \widetilde K_t+(1-\theta)\widetilde H_t-\widetilde Y_t\\
&0\approx \widetilde Y_t-\widetilde K_t-\widetilde r_t
\end{align*}

随机过程$\lambda_{t+1}=\gamma \lambda_t +\varepsilon_{t+1} $，也可以对数线性化后，为
\[
\widetilde \lambda_{t+1}=\gamma \widetilde \lambda_t +\mu_{t+1}
\]
其中，$\bar \lambda=1,\mu_{t+1}=\varepsilon_{t+1}-(1-\gamma) $。
对数线性化后的方程组作为一个线性模型，就可以用矩阵来表达，如果把t期所有内生变量定义成一个向量形式，
\[x_t=[\widetilde K_{t+1}\quad \widetilde Y_t\quad \widetilde C_t\quad \widetilde H_t\quad \widetilde r_t]'\]
之所以是$\widetilde K_{t+1} $而不是$\widetilde K_t $，是因为$\widetilde K_t $ 是一个参数，对于$t$期而言是外生的而不是内生的变量了，这可以从值函数的表达上看出来。另一个外生的随机变量，
\[z_t=\widetilde \lambda_t\]
这样，Hansen模型的五个方程就可以用矩阵表达，
\[0=E_t[Fx_{t+1}+Gx_t+Hx_{t-1}+Lz_{t+1}+Mz_t]\]
其中，
\[F=
\begin{bmatrix}
0 & 0 & -1 & 0 & \beta \bar r\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
\end{bmatrix} ,
G=
\begin{bmatrix}
0 & 0 & 1 & 0 & 0\\
0 & 1 & -1 & -\frac{1}{1-\bar H} & 0\\
-\bar K & \bar Y & -\bar C & 0 & 0\\
0 & -1 & 0 & 1-\theta & 0\\
0 & 1 & 0 & 0 & -1\\
\end{bmatrix}
\]
\[
H=
\begin{bmatrix}
0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0\\
\bar K(1-\delta) & 0 & 0 & 0 & 0\\
\theta & 0 & 0 & 0 & 0\\
-1 & 0 & 0 & 0 & 0\\
\end{bmatrix},L=\begin{bmatrix}
0 \\ 0 \\ 0\\ 0 \\ 0
\end{bmatrix},
M=\begin{bmatrix}
	0 \\ 0 \\ 0 \\ 1 \\ 0
\end{bmatrix}\]

随机过程可以写为，
\[z_{t+1}=Nz_t+\mu_{t+1}\]
其中$E_t(\mu_{t+1})=0 $，$N=[\gamma] $。要求解上述模型，无非就是要得到如下均衡运动方程，
\[x_t=Px_{t-1}+Qz_t\]

Uglig(1999)定理一表明，如果该方程的解存在，那么矩阵$P$可以通过解如下矩阵二次方程而得到，
\[0=FP^2+GP+H\]

而$Q$由$Vvec(Q)=-vec(LN+M) $给出，其中$V=N'\otimes F+I_k\otimes(FP+G) $，这里，$I_k $是一个$k$ 维的单位向量，$k$的取值等于随机变量的个数（即$z_t $的维数）。

为什么呢？可将$x_t=Px_{t-1}+Qz_t $代入$0=E_t[Fx_{t+1}+Gx_t+Hx_{t-1}+Lz_{t+1}+Mz_t] $，然后整理可得，
\[0=(FP^2+GP+H)x_{t-1}+(FPQ+FQN+GQ+LN+M)z_t\]

由于该方程要求对于$x_{t-1} $和$z_t $的任何值都成立，所以括号内的项必须为0。解二次矩阵方程的根是复杂的。我们也许需要另外一种特殊的解法。

\subsection{运用跳跃变量的解法}
该解法的特殊性在于通过选择某一个变量为内生变量$x_t $，再将其他所有变量放入一个新的向量$y_t $中。这样做的好处是使得矩阵$P$成为了一个元素。

在上面例子中，可以选择
\[x_t=[\tilde K_{t+1}]\]

然后，
\[y_t=[\tilde Y_t\quad \tilde C_t\quad \tilde H_t\quad \tilde r_t]\]

然后将方程组按是否有期望项分组，模型的线性形式就可以写成，
\begin{align*}
& 0=Ax_t+Bx_{t-1}+Cy_t+Dz_t\\
& 0=E_t(Fx_{t+1}+Gx_t+Hx_{t-1}+Jy_{t+1}+Ky_t+Lz_{t+1}+Mz_t)\\
& z_{t+1}=Nz_t+\varepsilon_{t+1}\quad E_t(\varepsilon_{t+1}=0)
\end{align*}

具体到我们的模型，有，

\begin{align*}
& A=[0 \quad -\bar K \quad 0 \quad 0]',\\
& B=[0 \quad (1-\delta)\bar K \quad \theta \quad -1]'\\
& D=[0 \quad 0 \quad 1 \quad 0]'\\
& F=[0],G=[0],H=[0]\\
& J=[0 \quad -1 \quad 0 \quad \beta\bar r]\\
& K=[0 \quad 1 \quad 0 \quad 0 ]\\
& L=[0],M=[0],N=[\gamma]
\end{align*}

\[C=
\begin{bmatrix}
1 & -1 & -\frac{1}{1-\bar H} & 0 \\
\bar Y & -\bar C & 0 & 0\\
-1 & 0 & 1-\theta & 0 \\
1 & 0 & 0 & -1\\
\end{bmatrix}\]
该模型其解应为如下形式，
\begin{align*}
& x_t=Px_{t-1}+Qz_t\\
& y_t=Rx_{t-1}+Sz_t
\end{align*}

Uglig(1999)表明如果解存在(见书上的附录1有一个简单的推导过程)，必有，
\begin{align*}
0 &=(F-JC^{-1}A)P^2-(JC^{-1}B-G+KC^{-1}A)P-KC^{-1}B+H\\
R & =-C^{-1}(AP+B)
\end{align*}

同时$Q$和$S$满足，
\begin{align*}
 [N'\otimes(F-JC^{-1}A)+I_k\otimes(JR+FP+G-KC^{-1}A)]vec(Q)&=vec[(JC^{-1}D-L)N+KC^{-1}D-M]\\
 S&=-C^{-1}(AQ+D)
\end{align*}
其中$I_k $是单位阵，矩阵$Q$有$k$列。

\subsection{对数线性模型的校准}
有些参数，其经济含义十分明确，如$\theta $是国民收入中资本收入的占比，有些参数如$\beta,\delta $有一些被普遍接受的值，剩下的一些参数如$A,\gamma,\mu_{t+1} $等需要一些规则来确定。这些变量的选择应该使得模型中其他相关变量由此得到的推算值近似于实际数据中的观测值。譬如$A$的取值应该使得劳动时间占总时间的1/3，而$\gamma,\mu_{t+1} $应该使得$y_t $的方差和协方差与现实数据相匹配。

Hansen模型根据季度数据来校准，通过选择$\beta=0.99,\delta=0.025,\theta=0.36 $，而前面\eqref{eq9}式给出了$\overline H $的稳态，
\[\bar H=\frac{1}{1+\frac{A}{1-\theta}\left[1-\frac{\beta\delta\theta}{1-\beta(1-\delta)} \right]}\]

为了使得$ \overline{H} $为1/3，$A$的取值为2是合适的。实际上，我们取了$A=1.72$，此时$\overline H=0.3335 $。

为了获得$\gamma $，可从生产函数的实际数据估计得出，即
\[ln\lambda_t=lnY_t-\theta lnK_t-(1-\theta)H_t\]

得到$\{\lambda_t\} $序列一阶滞后自相关系数约为0.95，故$\gamma=0.95 $ 。

确定$\mu_{t+1} $的分布比较复杂，一般地，该分布方差应该使得模型模拟的产出方差与观察数据十分接近。要确定这一过程，只有得到了运动方程才能进行。

\eqref{eq10}式给出了稳态$\overline K $的表达，
\[\overline K=\overline H\left[ \frac{\theta \overline \lambda}{\frac{1}{\beta}-(1-\delta)}\right]\]

由于$A=1.72$，$\overline H=0.3335 $，从而$\overline K=12.6695 $，于是根据$\overline Y=\overline K^{\theta}\overline H^{1-\theta}=1.2353 $，而预算约束又表明$\bar C=\overline Y-\delta \overline K=0.9186 $，而$\frac{1}{\beta}=\overline r+(1-\delta)=\theta \frac{\overline Y}{\bar K}+(1-\delta) $，可得到$\overline r=0.0351 $。将这些值代入前面矩阵$ A\sim N$，可以得到一个关于$P$的二次方程，
\[0=7.0734\cdot P^2-14.2376\cdot P+7.1448\]

解之可知$P=0.9537$可以使系统稳定均衡。得到了$P$，就很容易求出$Q,R$和$S$。于是可以得到在稳态附近的五个运动方程，
\begin{align*}
\widetilde K_{t+1} & =0.9537\widetilde K_t+0.1132\widetilde \lambda_t\\
\widetilde Y_{t} & =0.2045\widetilde K_t+1.4523\widetilde \lambda_t\\
\widetilde C_{t} & =0.5691\widetilde K_t+0.3920\widetilde \lambda_t\\
\widetilde H_{t} & =-0.2430\widetilde K_t+0.7067\widetilde \lambda_t\\
\widetilde r_{t} & =-0.7955\widetilde K_t+1.4523\widetilde \lambda_t
\end{align*}
\subsection{模型中变量的方差}
数据中产出对数的标准差为1.76\%，而模型中的变量$\widetilde Y_t=lnY_t-ln\overline Y $，可知$\widetilde Y_t $的标准差也为1.76\%。

现在可以使用2个运动方程，从$ \widetilde Y_t $的方差出发来解决$ u_{t+1} $的概率分布问题。
\begin{align}
 \widetilde{K}_{t+1}=a\widetilde{K_t}+b\widetilde{\lambda_t} \label{abc_eq11}\\ 
 \widetilde{Y_{t}}=c\widetilde{K_t}+d\widetilde{\lambda_t} \label{abc_eq12}
\end{align}
再加一个技术的随机过程：
\[ \widetilde{\lambda_{t}}=\gamma\widetilde{\lambda}_{t-1}+\varepsilon_t  \]

具体运算如下：
\begin{enumerate}
	\item 将随机过程代入资本运动方程，
	\[ \widetilde{K}_{t+1}=a\widetilde{K}_t+b\gamma \widetilde{\lambda}_{t-1}+b\varepsilon_t \]
	\item 对上式进行无限递归代入，得，
	\[ \widetilde{K}_{t+1}=b\sum_{i=0}^{\infty}\sum_{j=0}^{i}a^j\gamma^{i-j}\varepsilon_{t-i} \quad\text{技术水平也能无限递归，}\widetilde{\lambda_t}=\sum_{i=0}^{\infty} \gamma^i\varepsilon_{t-i}\]
	\item 将上两式代入产出方程\eqref{abc_eq12}式，得，
	\[ \widetilde{Y}_t=(d\cdot \varepsilon_t+\sum_{i=0}^{\infty} [cb\sum_{j=0}^i a^j\gamma^{i-j}+d\gamma^{i+1}])\varepsilon_{t-1-i} \]
	因为技术冲击相互独立，那么方差为，
	\[ \widetilde{Y}_t=(d^2+\sum_{i=0}^{\infty} [cb\sum_{j=0}^i a^j\gamma^{i-j}+d\gamma^{i+1}]^2)Var(\varepsilon_t) \]
	括号中的项是收敛的，因此可以得到任何想要的精度。
	\item 当得到了$\varepsilon_t$的方差后，可以观察到\eqref{abc_eq12}式与消费、工作时间或者资本租金的运动方程是类似的，那么就可以通过使用他们的运动方程(即不同的$c,d$)与\eqref{abc_eq11}式类似地联合，来得到消费、工作时间或者资本租金的方差。
\end{enumerate}

书中的表1是由模型得到的各变量相对产出的方差，表2是真实世界中各变量的相对产出的方差，足见，该模型对消费的近似是合理的，但对投资和工作时间的近似是相当离谱的。

\section{含有不可分劳动的 Hansen 模型}
考虑$ t $期居民以$ \alpha_t $的概率向厂商提供$ h_0 $单位的劳动，从而厂商的劳动力需求为$ h_t=\alpha h_0 $，居民在$ t $期是否工作随机决定。于是，时期$ t $的期望效用为，
\[ u(c_t,\alpha_t)=\ln c_t+\frac{h_t}{h_0}\cdot A\ln (1-h_0)+(1-\frac{h_t}{h_0})\cdot A\ln 1 \]

因为$ \ln 1=0 $，再进一步定义$ B=\frac{A\ln(1-h_0)}{h_0} $,则上式可以进一步简化为，
\begin{align*}
 u(c_t,\alpha_t)&=\ln c_t+\frac{h_t}{h_0}\cdot A\ln (1-h_0)\\
 & = \ln c_t+Bh_t
\end{align*}
那么，最大化问题为，
\begin{align*}
\max&\sum_{t=0}^{\infty}\beta^t[\ln c_t+Bh_t]\\
s.t. \quad& c_t+i_t=w_th_t+r_tk_t\\
& k_{t+1}=(1-\delta)k_t+i_t\\
& \ln\lambda_{t+1}=\gamma\ln\lambda_t+\varepsilon_{t+1}
\end{align*}

从社会计划者的角度来看，其最大化问题为(解与上述竞争均衡解一致)，
\begin{align*}
\max&\sum_{t=0}^{\infty}\beta^t[\ln c_t+Bh_t]\\
s.t. \quad& \lambda_tk_t^{\theta}h_t^{1-\theta}=c_t+k_{t+1}-(1-\delta)k_t\\
& \ln\lambda_{t+1}=\gamma\ln\lambda_t+\varepsilon_{t+1}
\end{align*}


\textbf{以Bellman方程的形式写为}，
\begin{align*}
 V(k_t\lambda_t)=& \max_{c_t,h_t} [\ln c_t+Bh_t+\beta E_tV(k_{t+1},\lambda_{t+1})] \\
s.t. \quad& \lambda_tk_t^{\theta}h_t^{1-\theta}=c_t+k_{t+1}-(1-\delta)k_t\\
& \ln\lambda_{t+1}=\gamma\ln\lambda_t+\varepsilon_{t+1}
\end{align*}

为方便求解，将控制变量更换为$ k_{t+1},h_t $,并进行代入，可以得到Bellman方程的另一种表达，
\[ V(k_t\lambda_t)=\max_{c_t,h_t} [\ln (\lambda_tk_t^{\theta}h_t^{1-\theta}-k_{t+1}+(1-\delta)k_t)+Bh_t+\beta E_tV(k_{t+1},\lambda_{t+1})] \]

上式的一阶条件中$ \lambda_tk_t^{\theta}h_t^{1-\theta}-k_{t+1}+(1-\delta)k_t $部分用$ c_t $替代，则有，
\begin{align*}
0&=\frac{1}{c_t}[(1-\theta)\lambda_tk_t^{\theta}h_t^{-\theta}]+B\\
0&= -\frac{1}{c_t}+E_t\left[ \frac{1}{c_{t+1}}\theta\lambda_{t+1}^{\theta-1}h_t^{1-\theta}+(1-\delta) \right]
\end{align*}

由于家庭完全同质，则可用大写字母替换小写字母作为一种总量替换，于是一阶条件为，
\begin{align*}
C_t&=-\frac{(1-\theta)Y_t}{BH_t}\\
1&=\beta E_t\left[\frac{C_t}{C_{t+1}}(r_{t+1}+(1-\delta))\right]
\end{align*}

为封闭模型，再加上流量预算约束、生产函数和两个要素市场条件\footnote{关于工资的要素条件由于零利润条件使得可以用$ Y_t $替换$ w_tH_t+r_tK_t $，因此，关于工资的方程也成为冗余。}，
\begin{align*}
C_t+K_{t+1}&=Y_t+(1-\delta)K_t\\
r_t&=\theta \lambda_tK_t^{\theta-1}H_t^{1-\theta}\\
Y_t&=\lambda_tK_t^{\theta}H_t^{1-\theta}
\end{align*}

\subsection{稳态}
稳态即要求$ \overline{X}=X_t=X_{t+1}$,技术的稳态$\bar \lambda=1 $,于是关于一阶条件和各约束就简化为，
\begin{align*}
\frac{1}{\beta}&=\bar r+(1-\delta)\\
\overline{C}&=-\frac{(1-\theta)\overline{Y}}{B\overline{H}}\\
\bar r&=\theta \overline{K}^{\theta-1}\overline{H}^{1-\theta}\\
\overline{Y}&=\overline{K}^{\theta}\overline{H}^{1-\theta}\\
\overline{C}&=\overline{Y}-\delta\overline{K}
\end{align*}
\subsubsection{对数线性化}
对一阶条件和各约束对数线性化，可得到，
\begin{align*}
0&\approx \widetilde{C}_t-E_t\widetilde{C}_{t+1}+\beta \bar r E_t\widetilde{r}_{t+1}\\
0&\approx \widetilde{C}_t+\widetilde{H}_t-\widetilde{Y}_t\\
0&\approx \overline{Y}\widetilde{Y_t}-\overline{C}\widetilde{C}_t+(1-\delta\overline{K}\widetilde{K}_t)-\overline{K}\widetilde{K}_{t+1}\\
0&\approx \widetilde{Y}_t-\widetilde{\lambda}_t-\theta \widetilde{K}_t-(1-\theta)\widetilde{H}_t\\
0&\approx \widetilde{Y}_t-\widetilde{K}_t-\widetilde{r}_t
\end{align*}

在仅有一个条件包含期望的情况下，就可以将上述模型写成如下形式，
\begin{align*}
&0=Ax_t+Bx_{t+1}+Cy_t+Dz_t\\
&0=E_t(Fx_{t+1}+Gx_t+Hx_{t-1}+Jy_{t+1}+Ky_{t}+Lz_{t+1}+Mz_t)\\
&z_{t+1}=Nz_t+\varepsilon_{t+1}\qquad E_t(\varepsilon_{t+1})=0
\end{align*}

其中，各变量
\[x_t=[\widetilde{K}_{t+1}],\quad
y_t=\begin{bmatrix}
\widetilde{Y}_t\\\widetilde{C}_t\\\widetilde{H}_t\\\widetilde{r}_t
\end{bmatrix},\quad
z_t=[\widetilde{\lambda}_t]\]

各系数，
\[A=\begin{bmatrix}
0\\-\overline{K}\\0\\0
\end{bmatrix},
B=\begin{bmatrix}
0\\ \overline{K}(1-\delta)\\ \theta\\-1
\end{bmatrix},
C=\begin{bmatrix}
1 & -1 & -1&0\\
\overline{Y}& -\overline{C}& 0& 0\\
-1& 0& 1-\theta&0\\
1& 0& 0&-1
\end{bmatrix},
D=\begin{bmatrix}
0\\0\\1\\0
\end{bmatrix},
J=\begin{bmatrix}
0\\-1\\0\\ \beta \bar r
\end{bmatrix},
K=\begin{bmatrix}
0\\1\\0\\0
\end{bmatrix}\]
\[F=[0],G=[0],H=[0],L=[0],M=[0],N=[\gamma]\]

该问题的解是一个矩阵集$ P,Q,R,S $，
\begin{align*}
x_t=Px_{t-1}+Qz_t\\y_t=Rx_{t-1}+Sz_t
\end{align*}

与上一节类似的方法可以得到这些矩阵，从而，
\begin{align*}
\widetilde{K}_t=0.9418\widetilde{K}_{t-1}+0.1552\lambda_t\\
y_t=R\widetilde{K}_{t-1}+S\lambda_t
\end{align*}
其中，
\[R=
\begin{bmatrix}
0.0550\\0.5316\\-0.4766\\-0.9450
\end{bmatrix},
S=\begin{bmatrix}
1.9418\\0.4703\\1.4715\\1.9417
\end{bmatrix}
\]
\section{脉冲响应函数}
给模型的随机误差一个冲击时，模型会做出相应反应。脉冲响应的计算步骤如下：
\begin{enumerate}
	\item 在Hansen模型中，仅有技术包含随机冲击$ \varepsilon_t $，即
	\[\tilde \lambda_t=\gamma\tilde{\lambda}_{t-1}+\varepsilon_t  \]
	\item 当给$ \varepsilon_t $一个冲击，那么那个“单独的变量”遵循如下路径，
	\[ \tilde{K}_{t+1}=P\tilde{K}_t+Q\tilde{\lambda}_t\]
	\item 其他内生变量$ y_t=(\tilde{Y}_t,\tilde{C}_t,\tilde{H}_t,\tilde{r}_t) $遵循如下路径，
	\[ y_t=R\tilde{K}_t+S\tilde{\lambda}_t \]
\end{enumerate}
脉冲时从稳态开始，即意味着对数偏离从0开始。
\begin{figure}[H]
\includegraphics[scale=0.7]{irf.png}
\caption{脉冲响应图}
\small 源码见\lstinline|irf.R|。
\end{figure}

\section{Blanchard and Kahn的解法}
该解法是Dynare在一阶求解用到的解法，而且在极大似然估计时，也和该解法有关(李向阳(2018, p60)有一个关于线性解如何表达成状态空间形式的推导)。因此，这里再详细讲解该解法。一个线性随机差分方程可以写成状态空间形式，
\begin{equation}\label{hansen_eq_bc}
B\begin{bmatrix}
	x_{t+1}\\ E_ty_{t+1}
\end{bmatrix}=A\begin{bmatrix}
	x_t\\ y_t
\end{bmatrix}+G\varepsilon_t	
\end{equation}

其中$ x_t $是由前定变量构成的向量，$ y_t $是由非前定变量构成的向量，区别就是前定变量当期值(如$ t+1 $期)不由当期($ t+1 $期)扰动决定，但非前定变量当期值受当期扰动影响，这也是它为什么要用期望符号的原因。

\eqref{hansen_eq_bc}式可以重新写成，
\[ \begin{bmatrix}
	x_{t+1}\\ E_ty_{t+1}
\end{bmatrix}=B^{-1}A\begin{bmatrix}
	x_t\\ y_t
\end{bmatrix}+B^{-1}G\varepsilon_t	
 \]
 
 如果令$ Z=B^{-1}A $，那么可以对$ Z $进行特征值分解$ Z=\bar{M}\bar{\Lambda}\bar{M}^{-1} $，这里$ \bar{\Lambda} $是从小到大排序后的特征值，$ \bar{M} $是对应的特征向量。Blanchard and Kahn的解法要求的条件就是\textbf{大于1的特征值的数目要等于有期望符号变量的数目}，后面我们会看到这一点的重要性。
 
 \paragraph{只考虑确定性部分}先忽略随机项，那么模型重新写成，
 \begin{equation}\label{hansen_eq_xy}
 \begin{bmatrix}
 	x_{t+1}\\ E_ty_{t+1}
 \end{bmatrix}=\bar{M}\bar{\Lambda}\bar{M}^{-1}\begin{bmatrix}
 	x_t\\ y_t
 \end{bmatrix}\Longrightarrow \bar{M}^{-1}\begin{bmatrix}
 	x_{t+1}\\ E_ty_{t+1}
 \end{bmatrix}=\bar{\Lambda}\bar{M}^{-1}\begin{bmatrix}
 	x_t\\ y_t
 \end{bmatrix}	
 \end{equation}
 
 若令，
 \[ \bar{M}^{_1}=\begin{bmatrix}
 	\hat{M}_{11}&\hat{M}_{12}\\
 	\hat{M}_{21}&\hat{M}_{22}
 \end{bmatrix}, \bar{\Lambda}=\begin{bmatrix}
 \bar{\Lambda}_{11}&0\\
 0& \bar{\Lambda}_{22}
\end{bmatrix} \]
其中，$ \bar{\Lambda}_{11} $中包含的是模型中特征值小于1的，$ \bar{\Lambda}_{22} $是模型中特征值大于1的，即导致模型不稳定的特征值。这样，就可以把\eqref{hansen_eq_xy}式写成两个方程，
\begin{align*}
	[\hat{M}_{11}x_{t+1}+\hat{M}_{12}E_ty_{t+1}]=&\bar{\Lambda}_{11}[\hat{M}_{11}x_t+\hat{M}_{12}y_t]\\	[\hat{M}_{21}x_{t+1}+\hat{M}_{22}E_ty_{t+1}]=&\bar{\Lambda}_{22}[\hat{M}_{21}x_t+\hat{M}_{22}y_t]
\end{align*}

因为$ \bar{\Lambda}_{22} $中的特征值都大于1，而我们寻找的是稳定解，因此$ [\hat{M}_{21}x_t+\hat{M}_{22}y_t] $必须为0，否则模型发散。这就意味着，
\begin{equation}\label{hansen_eq_bk1}
y_t= -(\hat{M}_{22})^{-1}\hat{M}_{21}x_t	
\end{equation}

因此我们首先假设了模型是确定性的，因此，还有
\begin{equation}\label{hansen_eq_bk2}
E_ty_{t+1}= y_{t+1}= -(\hat{M}_{22})^{-1}\hat{M}_{21}x_{t+1} 	
\end{equation}
将\eqref{hansen_eq_bk1}和\eqref{hansen_eq_bk2}式带入\eqref{hansen_eq_xy}式，就有，
\[ x_{t+1}=[\hat{M}_{11}-\hat{M}_{12}(\hat{M}_{22})^{-1}\hat{M}_{21}]^{-1}\hat{\Lambda}_{11}[\hat{M}_{11}-\hat{M}_{12}(\hat{M}_{22})^{-1}\hat{M}_{21}]x_t \]

\paragraph{加入随机冲击}
遵循前面的步骤，模型重新写成，
 \begin{equation}\label{hansen_eq_sto}
\bar{M}^{-1}\begin{bmatrix}
		x_{t+1}\\ E_ty_{t+1}
	\end{bmatrix}=\bar{\Lambda}\bar{M}^{-1}\begin{bmatrix}
		x_t\\ y_t
	\end{bmatrix}	+\bar{M}^{-1}B^{-1}G[\varepsilon_t]
\end{equation}
或者写成分块的形式，


\chapter{线性二次动态规划}
我们想将目标函数近似成二次的，其一阶导数就是一个一次项，有利于动态的刻画。具体地，对于$ n\times 1 $的状态变量$ x_t $，$ m\times 1 $的控制变量$ y_t $，我们希望找到这样的贴现二次目标函数，
\[ \sum_{t=0}^{\infty} \beta^t (x'_tRx_t+y'_tQy_t+2y'_tWx_t) \]

受约束于，
\[ x_{t+1}=Ax_t+By_t \]

其中，$ R,A $均为$ n\times n $阶，$ Q $是$ m\times m $阶，$ W $是$ m\times n $阶，$ B $是$ n\times m $阶。这种建模方式要求：
\begin{itemize}
	\item 预算约束是线性的；
	\item 目标函数必须是二次的；
\end{itemize}

\section{Kydland and Prescott的解法}\label{chp7}
模型的一般形式可以写为，
\begin{align*}
\sum_{t=0}^{\infty}& \beta^tF(x_t,y_t)\\
s.t.\quad & x_{t+1}=G(x_t,y_t)=Ax_t+By_t
\end{align*}

假设预算约束$ G $为线性(如果非线性则将其代入目标函数，使得为线性)。

于是$ F(x_t,y_t) $的二阶泰勒展开为，
\begin{align}
F(x_t,y_t)\approx & F(\bar x,\bar y)+[F_x(\bar x,\bar y)'\quad F_y(\bar x,\bar y)']
\begin{bmatrix}
x_t-\bar x\\
y_t-\bar y
\end{bmatrix}\nonumber \\
& + [(x_t-\bar x)'\quad(y_t-\bar y)']
\begin{bmatrix}
\frac{F_{xx}(\bar x,\bar y)}{2} & \frac{F_{xy}(\bar x,\bar y)}{2}\\
\frac{F_{yx}(\bar x,\bar y)}{2} & \frac{F_{yy}(\bar x,\bar y)}{2}
\end{bmatrix}
\begin{bmatrix}
x_t-\bar x\\
y_t-\bar y
\end{bmatrix}\label{eq13}
\end{align}

为将上式写成二次项的形式，可以进一步改写符号，定义，
\[z_t=\begin{bmatrix}
	1\\x_t\\y_t
\end{bmatrix},\text{稳态时，}\bar z=\begin{bmatrix}
1\\ \bar x\\ \bar y
\end{bmatrix}\]

由于$ x_t,y_t $分别是$ k $维和$ l $维，于是$ z_t $就是$ 1+k+l $维，那么可以构造一个如下矩阵，
\[M=\begin{bmatrix}
	m_{11} & m_{12} & m_{13}\\
	m_{21} & m_{22} & m_{23}\\
	m_{31} & m_{32} & m_{33}
\end{bmatrix}\]

其中，$ m_{11} $是$ 1\times 1 $，$ m_{22} $是$ k\times k $，$ m_{33} $是$ l\times l $，其余矩阵在相应位置使得$ M $恰好为方阵。这样的设置就使得，
\[ z'_tMz_t=m_{11}+(m_{12}+m'_{21})x_t+(m_{13}+m'_{31})y_t+x'_tm_{22}x_t+x'_t(m_{23}+m'_{32})y_t+y'_tm_{33}y_t \]

将上式与\eqref{eq13}式比对，可知$ m_{11} $包含的是常数项，$ m_{12}=m'_{21},m_{13}=m'_{31} $包含的是一次项系数，$ m_{22},m_{33},m_{23}=m'_{32} $是二次项系数，$ M $是一个对称矩阵，
\begin{align*}
m_{11}=F(\bar x,\bar y)-\bar x'F_x(\bar x,\bar y)-\bar y'F_y(\bar x,\bar y)+\frac{\bar x'F_{xx}(\bar x,\bar y)\bar x}{2}+\bar x'F_{xy}(\bar x,\bar y)\bar y+\frac{\bar y'F_{yy}(\bar x,\bar y)\bar y}{2}
\end{align*}
\begin{align*}
m_{12}&=m'_{21}=\frac{F_x(\bar x,\bar y)'-\bar x'F_{xx}(\bar x,\bar y)-\bar y'F_{yx}(\bar x,\bar y)}{2}\\
	m_{13}&=m'_{31}=\frac{F_y(\bar x,\bar y)'-\bar x'F_{xy}(\bar x,\bar y)-\bar y'F_{yy}(\bar x,\bar y)}{2}
\end{align*}
\begin{align*}
m_{22}&= \frac{F_{xx}(\bar x,\bar y)}{2}\\
m_{23}&=m'_{32}= \frac{F_{xy}(\bar x,\bar y)}{2}\\
m_{33}&= \frac{F_{yy}(\bar x,\bar y)}{2}\\
\end{align*}

至此，我们就得到了一个一般化的二次贴现的动态规划问题，
\begin{align*}
&\sum_{t=0}^{\infty} \beta^tz'_tMz_t\\
s.t. \quad & \begin{bmatrix}
1\\x_{t+1}
\end{bmatrix}
=A\begin{bmatrix}
1\\ x_t
\end{bmatrix}
+By_t
\end{align*}

其中，$z'_t=[1\quad x_t\quad y_t]$。
\subsection{解Bellman方程}
上述目标函数的二次型可以具体为，
\begin{align}
z'_tMz_t&=[x'_t\quad y'_t]\begin{bmatrix}
R & W'\\
W & Q
\end{bmatrix}
\begin{bmatrix}
x_t\\y_t
\end{bmatrix}\nonumber\\
&= x'_tRx_t+y'_tQy_t+2y'_tWx_t\label{abc_eq14}
\end{align}

注意，此处$x_t $第一个元素为1。该模型的值函数可以表达为$ x'_tPx_t $，其中$ P $待定。如果$ P $存在，那么Bellman方程为，
\[ x'_tPx_t=\max_{y_t} (z'_tMz_t+\beta x'_{t+1}Px_{t+1}) \]

服从约束，
\[ x_{t+1}=Ax_t+By_t \]

由于$ x_t $第一个元素为1，因此，$A$第一行第一个的元素为1，该行其他元素为0。$ B $第一行全部元素为0。\textbf{将\eqref{abc_eq14}式代入Bellman方程，可得}，
\[ x'_tPx_t=\max_{y_t}[x'_tRx_t+y'_tQy_t+2y'_tWx_t+\beta(Ax_t+By_t)'P(Ax_t+By_t)] \]

其一阶条件为，
\begin{equation}\label{eq17}
(Q+\beta B'PB)y_t=-(W+\beta B'PA)x_t
\end{equation}


于是，政策函数为，
\[ y_t=Fx_t=-(Q+\beta B'PB)^{-1}(W+\beta B'PA)x_t \]

不过矩阵$ P $仍然待定，此时将政策函数代入Bellman方程以消掉$ y_t $，再经过大量的矩阵运算，就可以得到，
\[ P=R+\beta A'PA-(\beta A'PB+W')(Q+\beta B'PB)^{-1}(\beta B'PA+W) \]

然后通过给$ P $赋一个初值如$ P_0=I $，就可以不断迭代以得到收敛的序列$ \{P_k\} $，
\begin{equation}\label{eq15}
P_{k+1}=R+\beta A'P_kA-(\beta A'P_kB+W')(Q+\beta B'P_kB)^{-1}(\beta B'P_kA+W)
\end{equation}

\subsection{校准}
在参数$ \beta=0.99,\delta=0.025,\theta=0.36,A=1.72 $时，其稳态为$ \bar h=0.3335,\bar k=12.6695,\bar y=1.2353,\bar c=0.9186 $，于是，各矩阵如下，
\[a=
\begin{bmatrix}
-0.6056 & 0.5986 & -1.3823\\
0.5968 & -0.5926 & 1.4048\\
-1.3823 & 1.4048 & -6.6590
\end{bmatrix},M=\begin{bmatrix}
-1.6374 & 1.0996 & -1.0886 & 1.9361\\
1.0996 & -0.6056 & 0.5986 & -1.3823\\
-1.0886 & 0.5986 & -0.5926 & -1.4048\\
1.9361 & -1.3823 & 1.4048 & -6.6590
\end{bmatrix}  \]

由$ M=\left [\begin{smallmatrix}
R & W'\\W & Q
\end{smallmatrix}\right ] $，就可以得到矩阵$R,Q,W$，此时，就可以将单位阵$ P_0=\left[\begin{smallmatrix}
1 & 0\\0 &1
\end{smallmatrix}\right] $在式\eqref{eq15}中经历200次迭代就会逐渐稳定下来，得到，
\[ P_1=\begin{bmatrix}
-0.7515 & 0.9987\\0.9987 & -0.4545
\end{bmatrix},P_2=\begin{bmatrix}
-1.6909 & 0.8247\\ 0.8247 & -0.1924
\end{bmatrix},\cdots, P_{200}=\begin{bmatrix}
-96.3615 & 0.8779\\ 0.8779 & -0.0259
\end{bmatrix}\]


利用$ P $就可以得到政策函数，
\[ F=\begin{bmatrix}
0.5869 & 0.9537\\ 0.4146 & -0.0064
\end{bmatrix} \]

实际上，此时我们也可以验算一下$ F $算得对不对，因为知道$ \left[\begin{smallmatrix}
k_{t+1}\\h_t
\end{smallmatrix} \right]=F\left[\begin{smallmatrix}
1\\k_t
\end{smallmatrix} \right] $，该式子所有变量此时都是已知的，验证其是否大约相等即可。

\section{加入随机冲击}
此时，要注意：把非线性预算约束放在目标函数$F(x_t,y_t)$中，把随机项放入线性预算约束中，即冲击表现为，
\[ x_{t+1}=Ax_t+By_t+C\varepsilon_{t+1} \]

其中，$ \varepsilon_t $是一个独立同分布的随机变量，$E_t(\varepsilon_{t+1})=0,Var(\varepsilon)=\Sigma$，$ C $是$ m\times n $阶，$ m $是状态变量个数，$ n $是$ \varepsilon_{t+1} $的长度。这样，代表性行为人目标函数最大贴现期望值为，
\[ \max_{y_t}E_0\sum_{t=0}^{\infty}\beta^tF(x_t,y_t) \]

为得到均衡点附近的动态演化，就要要找到目标函数的二次近似，在稳态$ (\bar x,\bar y) $的邻域内为：
\[ F(x_t,y_t)=F(z_t)\approx z'_tMz_t=[1\quad x'_t]R\begin{bmatrix}
1\\x_t
\end{bmatrix}+y'_tQy_t+2y'_tW\begin{bmatrix}
1\\x_t
\end{bmatrix} \]

于是模型就可以重新书写为,
\begin{align*}
\max \quad & E_0\sum_{t=0}^{\infty}\beta^tz'_tMz_t\\
s.t.\quad& x_{t+1}=Ax_t+By_t+C\varepsilon_{t+1}
\end{align*}

此时，定义存在一个矩阵$ P $使得$ x'_tPx_t+c $是该模型的值函数，于是模型就可以写成递归形式(Bellman方程形式)，
\begin{align*}
x'_tPx_t+c & =\max_{y_t} [z'_tMz_t+\beta E_0(x'_{t+1}Px_{t+1}+c)]\\
s.t.\quad& x_{t+1}=Ax_t+By_t+C\varepsilon_{t+1}
\end{align*}

将预算约束代入Bellman方程期望项中，可得，
\begin{equation}\label{eq16}
 x'_tPx_t+c=\max_{y_t}[z'_tMz_t+\beta x'_tA'PAx_t+2x'_tA'PBy_t+\beta y'_tB'PBy_t+\beta E_0(\varepsilon'_{t+1}C'PC\varepsilon_{t+1})+\beta c] 
\end{equation}


此时考察上式中的期望项$ E_0(\varepsilon'_{t+1}C'PC\varepsilon_{t+1}) $，如果令方阵$ G=[g_{ik}]=C'PC $，那么就有，
\[E_t(\varepsilon'_{t+1}C'PC\varepsilon_{t+1})=\sum_j\sum_k E_t(\varepsilon_{t+1}^jg_{jk}\varepsilon_{t+1}^k)=\sum_j g_{jj}E_t(\varepsilon_{t+1}^j\varepsilon_{t+1}^j)\]

又对于$ k\ne j $时，恒有$ E_t(\varepsilon_{t+1}^k\varepsilon_{t+1}^j)=0 $,故有，
\[ E_t(\varepsilon'_{t+1}C'PC\varepsilon_{t+1}) = trace (C'PC\Sigma) \]

这是一个常数。将上式代入式\eqref{eq16}中，有，
\[ x'_tPx_t+c=\max_{y_t}[z'_tMz_t+\beta x'_tA'PAx_t+2x'_tA'PBy_t+\beta y'_tB'PBy_t+\beta trace(C'PC\Sigma)]+\beta c  \]

考虑到该方程两边常数项要与常数项相等，就可以得到两个结论，
\begin{align*}
c & =\frac{\beta trace(C'PC\Sigma)}{1-\beta}\\
x_t'Px'_t & =\max_{y_t}(z'_tMz_t+\beta x'_tA'PAx_t+2x'_tA'PBy_t+\beta y'_tB'PBy_t)\\
 & = \max_{y_t}(x'_tRx_t+y'_tQy_t+2y'_tWx_t+\beta x'_tA'PAx_t+2x'_tA'PBy_t+\beta y'_tB'PBy_t)
\end{align*}

该问题一阶条件为，
\[ (Q+\beta B'PB)y_t=-(W+\beta B'PA)x_t \]

可以看到其等价于确定性情况下的一阶条件式\eqref{eq17}。因此，政策函数也是一样的，
\[ y_t=Fx_t= -(Q+\beta B'PB)^{-1}(W+\beta B'PA)x_t\]

当$ k\rightarrow\infty $时，也可以得到一个收敛的$ P_{k+1} $,
\[P_{k+1}=R+\beta A'P_kA-(\beta A'P_kB+W')(Q+\beta B'P_kB)^{-1}(\beta B'P_kA+W)  \]

\textbf{区别于确定性模型，为了模拟随机冲击的时间路径，}理解下式即可，
\[ x_1=Ax_0+By_0+C\varepsilon_1 \]

即有了状态变量$ x_0 $，利用政策函数找到$ y_0 $，然后从相应分布中抽取一组随机变量$ \varepsilon_{1} $，于是就得到了下一期的状态变量，循环往复。

\subsection{一个例子}
以第\ref{sec1}章经济体为例，代表性行为人的模型为，
\begin{align*}
\max_{k_{t+1},h_t} & \sum_{0}^{\infty}\beta^t[\ln(\lambda_tk_t^{\theta}h_t^{1-\theta}+(1-\delta)k_t-k_{t+1})+A(1-h_t)]\\
s.t.\quad & k_{t+1}=k_{t+1}\\
& \lambda_{t+1}=(1-\gamma)+\gamma\lambda_t+\varepsilon_{t+1}
\end{align*}

这种随机项建模加入了常数以保证$ \varepsilon_t $的均值为0，方差有界。定义$ z_t=(x'_t\quad y'_t)' $，其中状态变量$ x_t=(1\quad k_t\quad \lambda_t)' $，控制变量为$ y_t=(k_{t+1}\quad h_t)' $。

在这样的定义下，预算约束$ x_{t+1}=Ax_t+By_t+C\varepsilon_{t+1} $就可以具体地写成，
\[\begin{bmatrix}
	1\\ k_{t+1}\\\lambda_{t+1}
\end{bmatrix}=
\begin{bmatrix}
1 & 0 & 0\\0 & 0& 0\\1-\gamma & 0 & \gamma
\end{bmatrix}\begin{bmatrix}
1\\k_t\\ \lambda_t
\end{bmatrix}+
\begin{bmatrix}
0 & 0\\1 & 0\\0 & 0
\end{bmatrix}\begin{bmatrix}
k_{t+1}\\h_t
\end{bmatrix}+\begin{bmatrix}
0\\ 0\\ 1
\end{bmatrix}\varepsilon_{t+1}\]

而目标函数(单期)的二阶泰勒展开为，
\begin{align*}
u(\cdot)\approx &  \ln(\bar \lambda\bar k^{\theta}\bar h^{1-\theta}-\delta\bar k)+A\ln(1-\bar h)\\
& + \frac{1}{\bar c}\left[\theta \frac{\bar y}{\bar k}+(1-\delta)\right](k_t-\bar k)+\frac{\bar y}{\bar c}(\lambda_t-\bar \lambda)-\frac{1}{\bar c}(k_{t+1}-\bar k)+\left[(1-\theta)\frac{1}{\bar c}\frac{\bar y}{\bar c}-\frac{A}{1-\bar h}\right](h_t-\bar h)\\
& +[k_t-\bar k\quad \lambda_t-\bar \lambda\quad k_{t+1}-\bar k\quad h_t-\bar h] \begin{bmatrix}
a_{11} & \hat a_{1\lambda} & a_{12} & a_{13}\\
\hat a_{\lambda 1} & \hat a_{\lambda\lambda} & \hat a_{\lambda 2} &\hat  a_{\lambda 3}\\
a_{21} & \hat a_{2\lambda} & a_{22} & a_{23}\\
a_{31} & \hat a_{3\lambda} & a_{32} & a_{33}
\end{bmatrix}
\begin{bmatrix}
k_t-\bar k\\ \lambda_t-\bar \lambda\\ k_{t+1}-\bar k\\ h_t-\bar h
\end{bmatrix}
\end{align*}

其中，$ a_{ij} $与\ref{chp7}节中的含义相同，下标中包含$ \lambda $表示 $ a_{ij} $是与$ \lambda_t $相关的。另外，此处与\ref{chp7}节中的不同之处在于增加了$ (\lambda_t-\bar \lambda) $项，各$ \hat a_{ij} $为，
\begin{align*}
\hat a_{1\lambda}& = \hat a_{\lambda 1}=-\frac{\bar y}{\bar c^2}\left[\theta\frac{\bar y}{\bar k}+1-\delta\right]+\theta\frac{\bar y}{\bar c\cdot \bar k}\\
\hat a_{\lambda\lambda}&=-\frac{\bar y^2}{\bar c^2}\\
\hat a_{2\lambda}& = \hat a_{\lambda2}= \frac{\bar y}{\bar c^2}\\
\hat a_{3\lambda} &=\hat a_{\lambda 3}=-(1-\theta)\frac{\bar y^2}{\bar c^2\bar h}+(1-\theta)\frac{\bar y}{\bar c\cdot\bar h}
\end{align*}

于是目标函数可以写为二次型，
\begin{align*}
\max_{y_t} & \sum_{t=0}^{\infty}z_t'Mz_t\\
s.t. & \quad x_{t+1}=Ax_t+By_t+C\varepsilon_{t+1}
\end{align*}
其中，$ M $是一个$ 5\times 5 $的矩阵，
\[M=\begin{bmatrix}
m_{11} & m_{12} & m_{13} & m_{14} & m_{15}\\
m_{21} & a_{11} & \hat a_{1\lambda} & a_{12} & a_{13}\\
m_{31} & \hat a_{\lambda 1} & \hat a_{\lambda\lambda} & \hat a_{\lambda 2} &\hat  a_{\lambda 3}\\
m_{41} & a_{21} & \hat a_{2\lambda} & a_{22} & a_{23}\\
m_{51} & a_{31} & \hat a_{3\lambda} & a_{32} & a_{33}
\end{bmatrix}\]

有了$ M $，就有了$ R,Q,W $，再根据约束条件就有了$ A,B,C $。于是，按照\ref{chp7}节同样的参数设置，再增加$ \gamma=0.95 $，根据方程，
\[ P_{k+1}=R+\beta A'P_kA-(\beta A'P_kB+W')(Q+\beta B'P_kB)^{-1}(\beta B'P_kA+W) \]

通过反复迭代，就能得到矩阵$ P $，
\[P= \begin{bmatrix}
-124.0532 & 1.0657 & 15.6762\\
1.0657 & -0.0259 & -0.1878\\
15.6762 & -0.1878 & -1.9963\\
\end{bmatrix} \]

然后使用$ y_t=Fx_t=-(Q+\beta B'PB)^{(-1)}(W+\beta B'PA)x_t $，可以得到政策函数，
\[ F=\begin{bmatrix}
-0.8470 & 0.9537 & 1.4340\\
0.1789 & -0.0064 & 0.2357\\
\end{bmatrix} \]

有了政策函数，就可以得到运动方程，
\[ \begin{bmatrix}
1\\k_{t+1}\\ \lambda_{t+1}
\end{bmatrix}=\begin{bmatrix}
1 & 0& 0\\0 & 0& 0\\ 0.05 & 0 & 0.95\\
\end{bmatrix}\begin{bmatrix}
1\\k_t\\ \lambda_t
\end{bmatrix}+\begin{bmatrix}
0 & 0\\1 & 0\\ 0& 0\\
\end{bmatrix}\begin{bmatrix}
-0.8470 & 0.9537 & 1.4340\\0.1789 & -0.0064 & 0.2357\\
\end{bmatrix}\begin{bmatrix}
1\\k_t\\ \lambda_t\\
\end{bmatrix}+\begin{bmatrix}
0\\ 0\\ 1\\
\end{bmatrix}\varepsilon_{t+1} \]

化简后，可得，
\[ \begin{bmatrix}
1\\k_{t+1}\\ \lambda_{t+1}
\end{bmatrix}=\begin{bmatrix}
1 & 0& 0\\-0.8470 & 0.9537& 1.4340\\ 0.05 & 0 & 0.95\\
\end{bmatrix}\begin{bmatrix}
1\\k_t\\ \lambda_t
\end{bmatrix}+\begin{bmatrix}
0\\ 0\\ 1\\
\end{bmatrix}\varepsilon_{t+1} \]

上式可以看成，
\[ x_{t+1}=\Psi x_{t}+C\varepsilon_{t+1} \]

那么将右边不断递归代入，就可以得到移动平均表达，
\[ x_{t+1}=\sum_{i=0}^{\infty}\Psi^iC\varepsilon_{t-i}+\Psi^{\infty}x_{-\infty} \]

其中$ \Psi^{\infty} $为$ i\rightarrow\infty $时的极限，那么有，
\[ \Psi^{\infty}=\begin{bmatrix}
1 & 0 &0\\ \bar k& 0 & 0\\\bar \lambda & 0 &0\\
\end{bmatrix}=\begin{bmatrix}
1 & 0&0\\12.6695 & 0& 0\\1 & 0& 0\\
\end{bmatrix} \]

注意到由于$ \Psi $的结构，使得任意$ x_{-\infty} $与$ \Psi $相乘都等于$ \bar x $。那么$ x $的协方差矩阵就可以写为，
\begin{align*}
 Var(x)&=E(x_{t+1}-\bar x)(x_{t+1}-\bar x)'  =E[(\sum_{i=0}^{\infty}\Psi^iC\varepsilon_{t-i})(\sum_{i=0}^{\infty}\varepsilon'_{t-i}C'\Psi'^i)]\\
 & = \sum_{i=0}^{\infty}\Psi^iCVar(\varepsilon_{t})C'(\Psi')^i\\
 & = Var(\varepsilon_{t})\sum_{i=0}^{\infty}\Psi^iCC'(\Psi')^i
\end{align*}

根据$ \Psi,C $的定义，就可以得到，
\[ Var(x) =Var\begin{bmatrix}
1\\ k_{t+1}\\ \lambda_{t+1}
\end{bmatrix} =Var(\varepsilon_t)\begin{bmatrix}
0 & 0& 0\\0& 4728.5 & 148.7\\0 & 148.7 & 10.3\\
\end{bmatrix}\] 

于是利用政策函数就可以得到控制变量的方差\footnote{原著似乎漏掉了$ Var(\varepsilon_t) $}，
\[ Var(y)=Var\begin{bmatrix}
k_{t+1}\\h_t
\end{bmatrix}=FVar(x)F'=Var(\varepsilon_t)\begin{bmatrix}
4728.2 & 6.7\\6.7 & 0.3\\
\end{bmatrix} \]

\chapter{CIA模型}
\section{基准框架}
在第\ref{inf}章做的那些非常一般化的推导是非常有用的。后续的模型都是在它们基础上的变化。如上一章的RBC模型，这儿的CIA模型(Cash in Advance)，就是对效用函数、生产函数的具体化。CIA更多地加入了货币。注意与第\ref{inf}章类似，大写字母表示总量。

每个家庭想最大化贴现后的预期效用函数，
\[ E_0\sum_{t=0}^\infty \beta^tu(c_t^i,h_t^i)=E_0\sum_{t=0}^\infty \beta^t\left\{ \ln c_t^i+\left[A\frac{\ln (1-h_0)}{h_0}\right]h_t^i\right\} \]

厂商有CD总量生产函数，
\[ y_t=\lambda_tK_t^\theta H_t^{1-\theta} \]
其中，
\[ \ln\lambda_{t+1}=\gamma\ln\lambda_t+\varepsilon_{t+1},\hspace{2em}\varepsilon_{t+1}\sim \mathcal{N}(0,\sigma_\varepsilon^2) \]

那么，从厂商的最优化来看，即工资等于劳动的边际生产力，利率等于资本的边际生产力，
\begin{align*}
 w_t=&(1-\theta)\lambda_tK_t^\theta H_t^{-\theta}\\
r_t=&\theta\lambda_tK_t^{\theta-1} H_t^{1-\theta}
\end{align*}
且劳动和资本的总量为，
\[ H_t=\int_{0}^1h_t^idi,\hspace{2em}K_t=\int_{0}^1k_t^idi \]

现在来看CIA约束。家庭将一定量的货币$ m_{t-1}^i $从上期持有至本期，并接受来自政府的转移支付$ (g_t-1)M_{t-1} $，其中$ M_{t-1} $是$ t-1 $期人均货币总量，因全部家庭的总量是1，所以人均货币总量就是总量。$ g_t $是$ t $期的货币增长率。因此，消费的CIA约束为，
\[ p_tc_t^i\le m_{t-1}^i+(g_t-1)M_{t-1} \]

要注意该模型满足条件$ g_t\ge\beta $。同时，家庭还面临流量约束，
\[ c_t^i+k_{t+1}^i+\frac{m_t^i}{p_t}=w_th_t^i+r_tk_t^i+(1-\delta)k_t^i+\frac{m_{t-1}^i+(g_t-1)M_{t-1}}{p_t} \]
该式左边是消费、新持有的资本存量和准备保留至下期的货币实际值存量。右边是工资、利息收入，折旧后的资本以及期初持有的货币量实际值。

当$ g=1 $时，自然是个稳态。当$ g\ne 1 $时，货币存量会随时间而变化，也就不存在$ p_t,m_t^i,M_t $的稳态。有两个办法可以解决该问题，一种是Cooley-Hansen的办法，将每个名义变量标准化，即定义，
\[ \hat{p}_t=\frac{p_t}{M_t},\hspace{2em}\hat{m}_t^i=\frac{m_t^i}{M_t},\hspace{2em}M_t/M_t=1 \]
从而CIA约束和流量约束可以写成，
\begin{align*}
&\frac{p_t}{M_t}c_t^i=\frac{m_{t-1}^i+(g_t-1)M_{t-1}}{M_t}=\frac{m_{t-1}^i+(g_t-1)M_{t-1}}{g_tM_{t-1}} \\
\Longrightarrow & \hat{p}_tc_t^i=\frac{\hat m_{t-1}^i+(g_t-1)}{g_t} \\
& c_t^i+k_{t+1}^i+\frac{\hat m_t^i}{\hat p_t}=w_th_t^i+r_tk_t^i+(1-\delta)k_t^i+\frac{\hat m_{t-1}^i+(g_t-1)}{g_t\hat p_t}
\end{align*}
这么写就可以让名义变量在$ g_t\ne 1 $时也拥有了稳态值。此时，我们重新归纳整个模型如下，
\begin{align*}
\max \hspace{1em} & E_0\sum_{t=0}^\infty \beta^t\left\{ \ln c_t^i+\left[A\frac{\ln (1-h_0)}{h_0}\right]h_t^i\right\}\\
& \text{两个流量约束条件：}\\
s.t.\hspace{1em}&c_t^i=\frac{\hat m_{t-1}^i+(g_t-1)}{ \hat{p}_tg_t} \\
& c_t^i+k_{t+1}^i+\frac{\hat m_t^i}{\hat p_t}=(1-\theta)\lambda_tK_t^\theta H_t^{-\theta}h_t^i+\theta\lambda_tK_t^{\theta-1} H_t^{1-\theta}k_t^i+(1-\delta)k_t^i+\frac{\hat m_{t-1}^i+(g_t-1)}{g_t\hat p_t}\\
&\text{两个外生冲击：}\\
&\ln\lambda_{t+1}=\gamma\ln\lambda_t+\varepsilon_{t+1}^\lambda\\
&\ln g_{t+1}=(1-\pi)\ln\bar{g}+\pi\ln g_t+\varepsilon_{t}^g,\;\;\text{这意味着}g_t\text{的稳态值为}\bar g\\
&\text{四个总量条件：}\\
&K_t=k_t^i,\;\;H_t=h_t^i\;\;C_t=c_t^i\;\; \hat{M}_t=\hat{m}_t^i=1
\end{align*}

该最优规划的拉格朗日函数为，
\begin{align*}
\mathcal{L} =\max_{c_t^i,h_t^i,k_{t+1}^i,\hat{m}_t^i} \hspace{1em}  & E_0\sum_{t=0}^\infty \beta^t \left[\ln c_t^i+Bh_t^i+\chi_t^1\left(\hat{p}_tc_t^i-\frac{\hat m_{t-1}^i+g_t-1}{ g_t}\right) \right. \\
+&\left. \chi_t^2\left(k_{t+1}^i+\frac{\hat{m}_t^i}{\hat{p}_t}-w_th_t^i-r_tk_t^i-(1-\delta)k_t^i\right) \right]
\end{align*}


其中，$ B=  A\frac{\ln (1-h_0)}{h_0}$。最优规划的一阶条件为，
\begin{align*}
\frac{\partial \mathcal{L}}{\partial c_t^i}&= \frac{1}{c_t^i}+\chi_t^1\hat{p}_t=0\\
\frac{\partial \mathcal{L}}{\partial h_t^i}&=B-\chi_t^2w_t=0\\
\frac{\partial \mathcal{L}}{\partial k_{t+1}^i}&= \chi_t^2-\beta E_t\chi_{t+1}[1-\delta+r_{t+1}]=0\\
\frac{\partial \mathcal{L}}{\partial \hat{m}_t^i}&=\frac{\chi_t^2}{\hat{p}_t}-\beta E_t\frac{\chi_{t+1}^1}{g_{t+1}}=0
\end{align*}
类似地，关键是先消去两个$ \chi $。这可以通过前两个一阶条件而搞定，$ \chi_t^1=-\frac{1}{\hat{p}_tc_t^i},\chi_t^2=\frac{B}{w_t}$，然后代入后两个一阶条件，再结合两个流量约束条件和要素定价公式，就有整个系统的方程组如下，
\begin{align*}
\frac{1}{\beta}&=E_t\frac{w_t}{w_{t+1}}(1-\delta+r_{t+1})\\
\frac{B}{w_t\hat{p}_t}&=-\beta E_t\frac{1}{\hat{p}_{t+1}c_{t+1}^ig_{t+1}}\\
\hat{p}_tc_t^i&=\frac{\hat m_{t-1}^i+(g_t-1)}{g_t}=1\\
k_{t+1}^i+\frac{\hat m_t^i}{\hat p_t}&=w_th_t^i+r_tk_t^i+(1-\delta)k_t^i\\
 w_t=&(1-\theta)\lambda_tK_t^\theta H_t^{-\theta}\\
r_t=&\theta\lambda_tK_t^{\theta-1} H_t^{1-\theta}
\end{align*}

注意模型的搭建，我们讲了一个经济学的故事，然后引入了三个全新的变量$ \hat{m}_t,M_t, \hat{p}_t $，在数学上要求解需要增加三个新的方程才能解决此问题。但你要注意到，在方程中，个体变量和总体变量存在确定的数量关系，这就意味着总可以把总体变量换算成个体变量，而在本文背景下，个体变量就等于总体变量。因此，我们实际上只引入了两个变量$ \hat{m}_t,\hat{p}_t $，但还要看到$ \hat{m}_t $和$ \hat{p} $在整个系统中总是可以以$ \frac{m_t}{p_t} $的形式成对出现，这意味着从比值识别出独立的$ m_t,p_t $在没有其他条件干预的情况下是困难的。但巧妙之处在于，因为个体值等于总体值（或者它们存在确定的比值关系也行），即$ \hat{M}_t=\hat{m}_t^i $，而$ \hat{M}_t=M_t/M_t=1 $，使得$ \hat{m}_t^i=1 $，从而相当于只引进了一个变量，当然，此时我们已经增加了一个关于$ \hat{m}_t^i $的一阶条件方程。因此，多了一个变量，多了一个方程，搞定。

然后来求解稳态。因为所有家庭相同，因此均衡时的个体值等于总量值，在稳态时有，
\[ \bar{C}=C_t=c_t^i,\bar{H}=H_t=h_t^i,\bar{K}=K_{t+1}=k_{t+1}^i,\bar{M}=\hat{M}=\hat{m}_t^i \]
于是上述六个方程构成的方程组在稳态下为，
\begin{align*}
\frac{1}{\beta}&=1-\delta+\bar{r}\\
\frac{B}{\bar{w}}&=-\frac{\beta}{\bar{C}\bar{g}}\\
\hat{p}\bar{C}&=1\\
\frac{1}{\hat p}&=\bar{w}\bar{H}+(\bar{r}-\delta)\bar{K}\\
\bar{w}=&(1-\theta)\left(\frac{\bar{K}}{\bar{H}}\right)^{-\theta}\\
\bar{r}=&\theta\left(\frac{\bar{K}}{\bar{H}}\right)^{\theta-1}
\end{align*}
\section{铸币税}
把政府以新发行货币购买的产品的实际价值定义为铸币税。在模型中，我们就不能简单地将货币分发给每个家庭，而是采取发行新货币的方式来填补财政赤字。这样的话，对前述基准模型做两个改动：
\begin{itemize}
	\item 一是要加入政府预算约束。即新发行的货币量要等于政府购买或者说政府赤字。政府的约束写成，
	\[ g_t=\hat{g}_t\bar{g}=\frac{M_t-M_{t-1}}{p_t} \]
	其中，$ g_t $是政府购买的实际水平，$ \hat{g} $是随机冲击，服从$ \ln \hat{g}_t=\pi\ln\hat{g}_{t-1}+\varepsilon_t^g ,\varepsilon_t^g\sim\mathcal{N}(0,\sigma_g^2)$，$\bar{g}$是政府购买或赤字的平均水平。可以看到左边是政府的实际购买量，右边就是新发行的货币量。如果定义$ \varphi_t $是$ t $期的货币增长率，那么有$ M_t=\varphi_tM_{t-1} $，然后与前面一样，对相关变量除以$ M_t $获得一个标准化的形式，即$ \hat{m}_t^i=m_t^i/M_t,\hat{p}=p_t/M_t $，那么政府预算可以重新写成，
	\[ g_t=\hat{g}_t\bar{g}=\frac{\frac{M_t}{M_t}-\frac{M_{t-1}}{M_t}}{\frac{p_t}{M_t}}=\frac{1-\frac{1}{\varphi_t}}{\hat{p}_t} \]
	\item 二是不再假定货币直接转移支付给家庭。
\end{itemize}


我们重新设定家庭的最优规划为，
\[\begin{aligned}
	\max_{c_t^i,h_t^i,k_{t+1}^i,\hat{m}_t^i}\;\; \sum_{t=0}^{\infty}\beta^tu(c_t^i,1-h_t^i)=\sum_{t=0}^{\infty}\beta^t(\ln c_t^i+Bh_t^i)\\
	s.t.\;\; \hat{p}_tc_t^i=\frac{\hat{m}_{t-1}^i}{\varphi_t}\\
	k_{t+1}^i+\frac{\hat{m}_t^i}{\hat{p}_t}=w_th_t^i+r_tk_t^i+(1-\delta)k_t^i
\end{aligned}
\]
其中第一个约束条件是家庭的现金约束，即当期的消费由上期保留至本期的货币支付\footnote{推到如下，\[ \begin{aligned}
		p_tc^i_t=m^i_{t-1} \\
		\Longrightarrow \frac{p_tc_t^i}{M_t}=\frac{m_{t-1}^i}{M_t}\\
		\Longrightarrow \hat{p}c_t^i=\frac{m_{t-1}^i/M_{t-1}}{M_t/M_{t-1}}=\frac{\hat{m}_{t-1}}{\varphi_t}
	\end{aligned}\]}。第二个约束条件是说新的资本存量和保留至下期的货币现金等于工资收入、利息收入和折旧后的资本。

针对该最优规划，可以把$ c_t^i,h_t^i $用约束条件替换掉，然后仅通过选择$ k_{t+1}^i,m_t^i $来优化目标函数，即
\[ \max_{\hat{m}_t^i,k_{t+1}^i} E_0\sum_{t=0}^{\infty}\beta^t\left[\ln\left(\frac{\hat{m}_{t-1}^i}{\varphi_t\hat{p}_t}\right)+B\left(\frac{k_{t+1}^i+\frac{\hat{m}_t^i}{\hat{p}}-r_tk_t^i-(1-\delta)k_t^i}{w_t}\right)\right]\]

对应的两个一阶条件为，
\begin{align*}
	\frac{1}{w_t}=\beta E_t\left(\frac{r_{t+1}+1-\delta}{w_{t+1}}\right)\\
	-\frac{B}{\hat{p}w_t}=\frac{\beta}{\hat{m}_t}
\end{align*}
对于整个系统而言，我们还需要$ w,r $的表达式。另外，如果你还感兴趣于$ c_t^i,h_t^i $，那么你要把这两个约束条件再放回来。同时还要给$ \varphi_t $一个行为方式，也就是政府的CIA约束。最后加上两个外生变量的随机运动，9个内生变量$ w_t,r_t,p_t,C_t,K_t,H_t,\varphi_t,\hat{g}_t,\lambda_t $构成的9个方程就可以构成整个系统如下，
\begin{align*}
	\frac{1}{w_t}=\beta E_t\left(\frac{r_{t+1}+1-\delta}{w_{t+1}}\right)\;\;\text{一阶条件}\\
	-\frac{B}{\hat{p}w_t}=\frac{\beta}{\hat{m}_t}=\beta\;\;\text{一阶条件}\\
	w_t=(1-\theta)\lambda_tK_t^\theta H_t^{-\theta}\;\;\text{工资}\\
	r_t=\theta\lambda_tK_t^{\theta-1}H_t^{1-\theta}\;\;\text{利率}\\
	\hat{p}_tC_t=\frac{\hat{m}_{t-1}^i}{\varphi_t}=\frac{1}{\varphi_t}\;\;\text{居民消费CIA}\\
	K_{t+1}+\frac{\hat{m}_t^i}{\hat{p}_t}=w_tH_t+r_tK_t+(1-\delta)K_t=K_{t+1}+\frac{1}{\hat{p}_t}\;\;\text{资源约束}\\
	\hat{p}_t\hat{g}_t\bar{g}=1-\frac{1}{\varphi_t}\;\;\text{政府CIA}\\
	\ln \hat{g}_t=\pi\ln\hat{g}_{t-1}+\varepsilon_t^g\\
	\ln\lambda_t=\gamma\ln\lambda_{t-1}+\varepsilon_t^\lambda
\end{align*}
针对该系统的稳态，可以获得$ \bar{g} $与$ \varphi $的关系，即铸币税与通胀率的关系。
\chapter{MIU}
我认为在数学形式上MIU比CIA更加简单些。MIU模型所做的就是把真实货币余额放入效用函数，其他没有变化。

家庭最优规划为，
\begin{align*}
	\max\;\;E_0\sum_{t=0}^{\infty}\beta^tu(c_t^i,\frac{m_t^i}{P_t},1-h_t^i)=\ln c_t^i+D\ln\frac{m_t^i}{P_t}+Bh_t^i\\
	s.t.\;\; c_t^i+k_{t+1}^i+\frac{m_t^i}{P_t}=w_th_t^i+r_tk_t^i+(1-\delta)k_t^i+\frac{m_{t-1}^i}{P_t}+(g_{t-1}-1)\frac{M_{t-1}}{P_t}
\end{align*}
其中，$ (g_{t-1}-1)\frac{M_{t-1}}{P_t} $是政府对居民的一次性货币转移。构造拉格朗日函数如下，
\[ \begin{aligned}
	L= &E_0\sum_{t=0}^\infty\beta^t\left[\ln c_t^i+D\ln\frac{m_t^i}{P_t}+Bh_t^i-\lambda_t\left(k_{t+1}^i-w_th_t^i-r_tk_t^i-(1-\delta)k_t^i-\frac{m_{t-1}^i}{P_t}-\right.\right.\\
	&\left.\left.(g_{t-1}-1)\frac{M_{t-1}}{P_t}+c_t^i+\frac{m_t^i}{P_t}\right)\right]\end{aligned}\]
针对$ c_t^i,m_t^i,h_t^i,k_{t+1}^i $的四个一阶条件如下，
\begin{align*}
	\frac{\partial L}{\partial c_t^i}=\frac{1}{c_t^i}-\lambda_t=0\\
	\frac{\partial L}{\partial m_t^i}=D\frac{P_t}{m_t^iP_t}-\frac{\lambda_t}{P_t}+\beta E_t\frac{\lambda_{t+1}}{P_{t+1}}=0\\
	\frac{\partial L}{\partial h_t^i}=B+\lambda_tw_t=0\\
	\frac{\partial L}{\partial k_{t+1}^i}=-\lambda_t+\beta E_t[\lambda_{t+1}(r_{t+1}+1-\delta)]=0
\end{align*}
通过消除拉格朗日乘子，有，
\begin{align}\label{eq_miu_m}
	\frac{1}{c_t^i}=\beta E_t\frac{P_t}{c_{t+1}^iP_{t+1}}+\frac{DP_t}{m_t^i}\\
	\frac{1}{c_t^i}=-\frac{B}{w_t}\\
	\frac{1}{c_t^i}=\beta E_t\frac{1}{c_{t+1}^i}[r_{t+1}+1-\delta]
\end{align}
再联合资源约束，$ w_t,r_t $的决定公式，
\begin{align}
	c_t^i+k_{t+1}^i+\frac{m_t^i}{P_t}=w_th_t^i+r_tk_t^i+(1-\delta)k_t^i+\frac{m_{t-1}^i}{P_t}+(g_{t-1}-1)\frac{M_{t-1}}{P_t}\\
	 w_t=(1-\theta)\lambda_tK_t^\theta H_t^{-\theta}\\\label{eq_miu_r}
	r_t=\theta\lambda_tK_t^{\theta-1} H_t^{1-\theta}
\end{align}
\eqref{eq_miu_m}-\eqref{eq_miu_r}式就构成了$ c_t^i,m_t^i,h_t^i,k_{t+1}^i,w_t,r_t $六个变量的运动系统。类似地，通过总量变量的替换，有，
\begin{align*}
	\frac{1}{C_t}=\beta E_t\frac{P_t}{C_{t+1}P_{t+1}}+\frac{DP_t}{M_t}\\
\frac{1}{C_t}=-\frac{B}{w_t}\\
\frac{1}{C_t}=\beta E_t\frac{1}{C_{t+1}}[r_{t+1}+1-\delta]\\
C_t+K_{t+1}+\frac{M_t}{P_t}=w_tH_t+r_tK_t+(1-\delta)K_t+\frac{M_{t-1}}{P_t}+(g_{t-1}-1)\frac{M_{t-1}}{P_t}\\
%\text{或者写成，}C_t+K_{t+1}+\overline{M/P}_t=w_tH_t+r_tK_t+(1-\delta)K_t+\frac{\overline{M/P}_{t-1}}{g_{t}}+(g_{t-1}-1)\frac{\overline{M/P}_{t-1}}{g_{t}}\\
w_t=(1-\theta)\lambda_tK_t^\theta H_t^{-\theta}\\
r_t=\theta\lambda_tK_t^{\theta-1} H_t^{1-\theta}\\
Y_t=\lambda_tK_t^\theta H_t^{1-\theta}\\
M_t=g_tM_{t-1}\\
\ln g_t=(1-\pi)\ln\bar{g}+\pi\ln g_{t-1}+\varepsilon_t^g\\
\ln\lambda_{t}=\gamma\ln\lambda_{t-1}+\varepsilon_{t}^\lambda
\end{align*}

注意以下几点，
\begin{itemize}
	\item 注意到方程中还多了个$ P_t $，需要一个新的方程，我们引入了$ M_t=g_tM_{t-1} $。
	\item 我们在最后增加了一个总产出的方程，对应地多了一个内生变量$ Y_t $，这主要可以在写代码计算产出脉冲响应时比较方便。	
	\item 在稳态时，$ M_t/P_t $会是一个常数，价格的增速即通胀率$ \bar{\pi} $与货币余额增速$ \bar{g} $相等。
\end{itemize}


\chapter{交错定价}
前面的模型对冲击的反应都过于迅速，回归稳态也过于迅速，对于价格更是如此。因此，我们需要迟滞价格。这就是Calvo的模型，每期随机选定部分企业可以调整价格。这里就引入了生产中间产品的企业。
\section{最终产品企业}
中间产品$ Y_t(k) $是一个连续统，用$ k\in [0,1] $表示，这些产品各不相同，但最终被某个完全竞争的产品企业打包生产成第$ t $期的最终产品$ Y_t $，这个生产函数可以写为，
\[ Y_t=\left[\int_{0}^1Y_t(k)^{\frac{\psi-1}{\psi}}dk\right]^{\frac{\psi}{\psi-1}} \]
其中$ \psi>1 $。很容易看到这是常替代弹性的生产技术。一个要利润最大化的厂商要解决的最优规划是，
\begin{equation}\label{eq_stg_py}
\max_{Y_t(k)}P_t\left[\int_{0}^1Y_t(k)^{\frac{\psi-1}{\psi}}dk\right]^{\frac{\psi}{\psi-1}}-\int_{0}^1P_t(k)Y_t(k)dk	
\end{equation}


其一阶条件为(可以把积分写成求和的形式，然后对于单个的$Y_t(k)$求导来理解)，
\[ P_t\left[\int_{0}^1Y_t(k)^{\frac{\psi-1}{\psi}}dk\right]^{\frac{1}{\psi-1}}Y_t(k)^{\frac{-1}{\psi}}=P_t(k) \]
这就简化成要素需求函数(这里体现为中间产品的需求)\footnote{\begin{align*}
	P_t\left[\int_{0}^1Y_t(k)^{\frac{\psi-1}{\psi}}dk\right]^{\frac{1}{\psi-1}}Y_t(k)^{\frac{-1}{\psi}}&=P_tY_t^{\frac{1}{\psi}}Y_t(k)^{\frac{-1}{\psi}}=P_t(k)\\
	\Longrightarrow Y_t(k)^{\frac{-1}{\psi}}& = \frac{P_t(k)}{P_tY_t^{\frac{1}{\psi}}}\\
	\Longrightarrow Y_t(k)&=Y_t\left(\frac{P_t}{P_t(k)}\right)^\psi
	\end{align*}}，
\[ Y_t(k)=Y_t\left(\frac{P_t}{P_t(k)}\right)^\psi \]

通过将这个单个中间产品的需求函数代入总生产函数，就可以得到价格指数或者说最终产品的定价规则，
\begin{align*}
& Y_t=\left[\int_{0}^1\left[Y_t\left(\frac{P_t}{P_t(k)}\right)^\psi\right]^{\frac{\psi-1}{\psi}}dk\right]^{\frac{\psi}{\psi-1}}=Y_t\left[\int_{0}^1\left(\frac{P_t}{P_t(k)}\right)^{\psi-1}dk\right]^{\frac{\psi}{\psi-1}}=Y_tP_t^\psi\left[\int_{0}^1\left(\frac{1}{P_t(k)}\right)^{\psi-1}dk\right]^{\frac{\psi}{\psi-1}}\\
 \Longrightarrow & P_t=\left[\int_{0}^1P_t(k)^{1-\psi}dk\right]^{\frac{1}{1-\psi}}
\end{align*}
\section{中间产品企业}
逻辑在于中间产品企业们不能每期选择自己产品的价格，只有$ 1-\rho $比例的企业随机选中，从而可以选择第$ t $期的价格，否则只能保持上一期价格不变。这就意味着，在$ t $期一旦价格确定，该价格就有$ \rho $的概率持续到$ t+1 $期，$ \rho^2 $的概率持续到$ t+2 $期。因此，具备在$ t $期选择价格的中间产品企业$ k $，它将如下最大化，
\begin{align}\nonumber
\max_{P_t(k)}&E_t\sum_{i=0}^\infty \beta^i\rho^i\left[P_t(k)Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi-P_{t+i}r_{t+i}K_{t+i}(k)-P_{t+i}w_{t+i}H_{t+i}(k)\right]\\\label{med1}
s.t.\;\;& Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi=\lambda_{t+i}K_{t+i}^\theta(k)H_{t+i}^{1-\theta}(k)
\end{align}

之所以在利率和工资前面要乘以价格$ P_{t+i} $是因为该模型我们是要引入货币因素的。在有货币因素的模型中，凡是涉及到价格的，名义变量都要乘以$ P_t $。但上述最优化并不意味着成本最小化。现在来看看成本最小化问题，
\begin{align*}
\min_{K_{t+i}(k),H_{t+i}(k)}&r_{t+i}K_{t+i}(k)+w_{t+i}H_{t+i}(k)\\
s.t.\;\;&Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi=\lambda_{t+i}K_{t+i}^\theta(k)H_{t+i}^{1-\theta}(k)
\end{align*}

类似地，从这两个一阶条件中消去拉格朗日乘子，有，
\[ \frac{(1-\theta)r_{t+i}}{\theta w_{t+i}}=\frac{H_{t+i}(k)}{K_{t+i}(k)} \]
利用上式和约束条件可以得到劳动和资本的需求为(注意，这里的$ \lambda $是技术冲击不是拉格朗日乘子)，
\begin{align}\label{eq_stg_h}
H_{t+i}(k)&= \frac{Y_{t+i}(k)}{\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta\\\label{eq_stg_k}
K_{t+i}(k)&= \frac{Y_{t+i}(k)}{\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^{\theta-1}
\end{align}

上面二式就是条件要素需求函数，将其带入成本最小化的目标函数就能得到最小化的成本为，
\begin{align}
& r_{t+i}\frac{Y_{t+i}(k)}{\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^{\theta-1}+w_{t+i}\frac{Y_{t+i}(k)}{\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta\\\label{eq_stg_cost}
\Longrightarrow & \frac{w_{t+i}}{(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta Y_{t+i}(k)
\end{align}

把这个最小化的成本代入到\eqref{med1}式的最优规划中，则有，
\begin{align*}
&\max_{P_t(k)}E_t\sum_{i=0}^\infty \beta^i\rho^i\left[P_t(k)Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi-\frac{P_{t+i}w_{t+i}}{(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta Y_{t+i}\left(\frac{P_{t+i}}{P_{t}(k)}\right)^\psi\right]\\
\Longrightarrow & \max_{P_t(k)}E_t\sum_{i=0}^\infty (\beta\rho)^iY_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi\left[P_t(k)-\frac{P_{t+i}w_{t+i}}{(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta\right]
\end{align*}

其一阶条件为，
\[ 0 = E_t\sum_{i=0}^\infty (\beta\rho)^iY_{t+i}(k)\left[1-\psi+\frac{\psi P_{t+i}w_{t+i}}{P_t(k)(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta\right]
 \]
从而得到中间产品的定价规则，
\[ P_t(k)=\frac{\psi}{1-\psi}\frac{E_t\sum_{i=0}^\infty (\beta\rho)^iP_{t+i}Y_{t+i}(k)\frac{w_{t+i}}{(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta}{E_t\sum_{i=0}^\infty (\beta\rho)^iY_{t+i}(k)} \]

\textbf{整个过程，我们再捋一遍，第一步，通过最终产品的利润最大化得到中间产品的需求函数，代入目标函数从而得到总价格指数。第二步，利用这个需求函数作为约束，通过中间产品的成本最小化得到劳动和资本的需求函数，代入目标函数从而得到成本函数。第三步，通过中间产品的利润最大化得到中间产品的定价公式。}\footnote{注意到，一个最优规划的求解次序，一般是得到一阶条件再消去拉格朗日乘子，然后一般还会代入原来的目标函数。这些步骤一般会产生具有相当经济内涵的表达。具体的，
\begin{itemize}
	\item 对于消费者的效用最大
	化问题(选择消费的数量以最大化效用)而言，这一步往往就得到了产品的需求函数，将这个需求函数代入效用函数，就得到了间接效用函数。
	\item 对于生产者利润最大化问题(选择生产的数量以最大化利润)而言，这就得到了劳动和资本的要素需求函数，代入目标函数就得到了产品的供给函数。
	\item 对于成本最小化问题(选择劳动或资本的数量以最小化成本)而言，这就得到了劳动和资本的条件要素需求函数，代入目标函数，就是企业的成本函数。
	\item 有时利润最大化问题和成本最小化问题是一个问题，有时不是。
\end{itemize}
}

\begin{figure}[H]
	\begin{tikzpicture}[inner sep=1mm,
	%	place/.style ={rectangle,fill =black!20, draw = black, minimum size = 1cm},
	place/.style ={rectangle, draw = black, minimum size = 1cm},
	arrow1/.style = {draw = black!50,-{Latex[length = 2mm, width = 1mm]}},]
	\node[place, label=above:A.最终产品利润最大化] (A) { $\max_{Y_t(k)}P_t\left[\int_{0}^1Y_t(k)^{\frac{\psi-1}{\psi}}dk\right]^{\frac{\psi}{\psi-1}}-\int_{0}^1P_t(k)Y_t(k)dk	$};
	\node[place] (B) [right = 2cm of A, label=above:B.得到中间产品需求函数]{$ Y_t(k)=Y_t\left(\frac{P_t}{P_t(k)}\right)^\psi$};
	\node[place] (F) [below = 1.5cm of A,  label = above:F.得到中间产品定价]{$P_t(k)=\frac{\psi}{1-\psi}\frac{E_t\sum_{i=0}^\infty (\beta\rho)^iP_{t+i}Y_{t+i}(k)\frac{w_{t+i}}{(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta}{E_t\sum_{i=0}^\infty (\beta\rho)^iY_{t+i}(k)}$};
	\node[place,align=center] (E) [below = 1.5cm of F,xshift = 2.5cm, label = above:E.中间产品利润最大化]{$ \max_{P_t(k)}E_t\sum_{i=0}^\infty \beta^i\rho^i\left[P_t(k)Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi-P_{t+i}r_{t+i}K_{t+i}(k)-P_{t+i}w_{t+i}H_{t+i}(k)\right]$\\
		$	s.t.\;\; Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi=\lambda_{t+i}K_{t+i}^\theta(k)H_{t+i}^{1-\theta}(k) $};
	\node[place, align=center] (C) [below = 1.5cm of E, xshift = -2.5cm, label=above:C.中间产品成本最小化]{$ \min_{K_{t+i}(k),H_{t+i}(k)} r_{t+i}K_{t+i}(k)+w_{t+i}H_{t+i}(k)$\\   $s.t.\;\;Y_{t+i}\left(\frac{P_{t+i}}{P_t(k)}\right)^\psi=\lambda_{t+i}K_{t+i}^\theta(k)H_{t+i}^{1-\theta}(k) $};
	\node[place, align=center] (D) [right = 1.5cm of C, label=above:D.得到劳动和资本的需求函数]{$ H_{t+i}(k)= \frac{Y_{t+i}(k)}{\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta$\\
	$	K_{t+i}(k)= \frac{Y_{t+i}(k)}{\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^{\theta-1} $};
	\node[place, fill = black!20] (G) [below = 1.62cm of B]{获得以资本和劳动表示的$ w $和$ r $};
	
\draw[arrow1] (A.east) to (B.west);
\draw[arrow1] (C.east) to (D.west);
\draw[arrow1] (E.west) to [bend left = 45] (F.west);
\draw[arrow1] (D.east) to [bend right = 45] node[right] {代入}(E.east);
\draw[arrow1, draw = blue, thick] (F.east) to (G.west);
\draw[arrow1, draw = blue, thick] (B.south) to (G.north);
\end{tikzpicture}
\caption{算法示意图}\label{fig-algo}
\end{figure}

现在，根据总的价格水平的定价规则，可以把价格的演变换种表达，
\[ P_t^{1-\psi}=\rho P_{t-1}^{1-\psi}+(1-\rho)P_t(k)^{1-\psi} \]

$ 1-\psi $次方是为了把价格变成水平形式，然后每期有$ \rho $比例不能调整价格，有$ 1-\rho $比例可以调整价格。
\section{家庭}
与CIA模型类似，家庭最优化问题为，
\begin{align*}
\max_{c_t,h_t}\quad &E_0\sum_{t=0}^\infty\beta^t(\ln c_t^i+Bh_t^i)\\
s.t.\quad & P_tc_t^i=m_{t-1}^i+(g_t-1)M_{t-1}\quad\text{CIA约束}\\
&k_{t+1}^i+\frac{m_t^i}{P_t}=w_th_t^i+r_tk_t^i+\xi_t^i+(1-\delta)k_t^i\quad\text{流量约束}
\end{align*}
和前面CIA模型唯一的区别在于这里多了个中间产品企业支付给家庭的超额利润$ \xi_t^i $。其两个一阶条件为，
\begin{align*}
	\frac{B}{w_t}=E_t\left[\frac{B\beta}{w_{t+1}}(r_{t+1}+1-\delta)\right]\\
	-E_t\left[\frac{\beta}{c_{t+1}^iP_{t+1}}\right]=\frac{B}{w_tP_t}
\end{align*}
接下来看，看看$ \xi_t^i $如何得到。根据\eqref{eq_stg_cost}式的总成本，可知中间产品的总的超额利润为，
\[ P_t\xi_t=\int_{0}^1P_t(k)Y_t(k)dk-P_t\frac{w_t}{(1-\theta)\lambda_t}\left[\frac{r_t(1-\theta)}{w_t\theta}\right]^\theta\int_{0}^1Y_t(k)dk \]
再依据\eqref{eq_stg_py}式最终企业的最优规划，因为它是完全竞争的，所以是零利润，这意味着，
\[ P_tY_t=\int_{0}^1P_t(k)Y_t(k)dk \]
综合上述两式，有加总后的超额利润，
\[ \xi_t= Y_t-\frac{w_t}{(1-\theta)\lambda_t}\left[\frac{r_t(1-\theta)}{w_t\theta}\right]^\theta\int_{0}^1Y_t(k)dk\]

\section{模型的全部方程}

一般而言，消费者问题可以得到关于$ C_t,K_t $的一阶条件，但往往要以$ w_t,r_t $作为已知条件。而$ w_t,r_t $是在企业模块的一阶条件中获得。那么对于非完全竞争情况，获得$ w_t,r_t $的表达更为复杂一些。在那里，往往是把最终产品生产设成完全竞争部门，中间产品生产设成垄断竞争部门。
\begin{itemize}
	\item 在最终产品生产中会得到中间产品的需求函数，即$ Y_t(k) $和$ P_t(k) $的函数关系，注意尽管该函数关系可能包含价格指数$ P_t $和加总产出$ Y_t $，但要意识到$ P_t $是$ P_t(k) $的函数，$ Y_t $是$ Y_t(k) $的函数。
	\item 然后，在中间产品生产中，因为是垄断竞争，这样的问题两步解决，第一步成本最小化，会得到$ w_t,r_t $关于$ Y_t(k) $的函数，第二步利润最大化，会得到，$ P_t(k) $关于$ Y_t(k) $的函数。那么这个函数和中间产品需求函数就可以解出$ P_t(k),Y_t(k) $，然后解出$ w_t,r_t $。	
\end{itemize}

从图\ref{fig-algo}来看，模块D给出$ w_t,r_t $关于$ Y_t(k) $的函数，模块F和B联立获得$ P_t(k) $和$ Y_t(k) $的表达，从而得到$ w_t,r_t $，然后联合消费者问题最终解出方程。



在这里总结模型方程如下，注意因为家庭同质，且为1个单位，因此，这里使用了总量条件$ C_t=c_t^i,K_t=k_t^i,H_t=h_t^i,M_t=m_t^i $。现在我们来捋一捋整个系统的方程构成。首先是消费者问题，消费者通过选择消费$ C_t $和$ H_t $来最优化，然后其有两个约束条件，因此消费者问题包含了四个方程，
\begin{align*}
	\frac{B}{w_t}=&E_t\left[\frac{B\beta}{w_{t+1}}(r_{t+1}+1-\delta)\right]\\
-E_t\left[\frac{\beta}{C_{t+1}P_{t+1}}\right]=&\frac{B}{w_tP_t}\\
P_tC_t=&M_{t-1}+(g_t-1)M_{t-1}\\
K_{t+1}+\frac{M_t}{P_t}=&w_tH_t+r_tK_t+\xi_t+(1-\delta)K_t
\end{align*}
注意到，这四个方程却包含了8个内生变量，即$ C_t,H_t,K_t,w_t,r_t,P_t,M_t,\xi_t $。因为厂商理论，根据我们前面几章的学习，知道它就是用来搞定$ w_t,r_t $的。在这里，厂商理论可以给出成本最小化的一阶条件，
\[ \frac{(1-\theta)r_t}{\theta w_t}=\frac{H_t(k)}{K_t(k)} \]
虽然这个一阶条件引入了新的变量$ H_t(k),K_t(k) $，但根据\eqref{eq_stg_h}式和\eqref{eq_stg_k}式，这两个新的变量都是$ Y_t(k),w_t,r_t $的函数，而，
\begin{align*}
Y_t(k)=\lambda_tK_t^\theta(k)H_t^{1-\theta}(k)\\
\ln\lambda_t=\gamma\ln\lambda_{t-1}+\varepsilon_t^\lambda
\end{align*}

因此，相当于模型前面5个内生变量$ C_t,H_t,K_t,w_t,r_t$都有了着落。接下来看后面3个内生变量，对于$ P_t $，可用如下两式搞定它，
\begin{align*}
 P_t^{1-\psi}=&\rho P_{t-1}^{1-\psi}+(1-\rho)P_t(k)^{1-\psi}\\
P_t(k)=&\frac{\psi}{1-\psi}\frac{E_t\sum_{i=0}^\infty (\beta\rho)^iP_{t+i}Y_{t+i}(k)\frac{w_{t+i}}{(1-\theta)\lambda_{t+i}}\left[\frac{r_{t+i}(1-\theta)}{w_{t+i}\theta}\right]^\theta}{E_t\sum_{i=0}^\infty (\beta\rho)^iY_{t+i}(k)}
\end{align*}

对于$ M_t $，它有如下方程，
\begin{align*}
M_t=&g_tM_{t-1}	\\
\ln g_t=&\pi\ln g_{t-1}+\varepsilon_t^g
\end{align*}
对于$ \xi_t $，有，
\begin{align*}
\xi_t=& Y_t-\frac{w_t}{(1-\theta)\lambda_t}\left[\frac{r_t(1-\theta)}{w_t\theta}\right]^\theta\int_{0}^1Y_t(k)dk\\
Y_t(k)=&Y_t\left(\frac{P_t}{P_t(k)}\right)^\psi\\
\int_{0}^1Y_t(k)dk=&\lambda_tK_t^\theta H_t^{1-\theta}
\end{align*}
最后一个等式可以从劳动力市场加总得到。因为，
\begin{align*}
	H_t=\int_{0}^1H_t(k)dk=\frac{1}{\lambda_t}\left[\frac{r_t(1-\theta)}{w_t\theta}\right]^\theta\int_0^1Y_t(k)dk\\		K_t=\int_{0}^1H_t(k)dk=\frac{1}{\lambda_t}\left[\frac{r_t(1-\theta)}{w_t\theta}\right]^{\theta-1}\int_0^1Y_t(k)dk\\
\end{align*}

\chapter{模拟矩估计}
\begin{note}
	文献来源：Rahmandad, H., R. Oliva and N.D. Osgood, eds. Analytical methods for dynamics modelers. 2015, The MIT Press: London.
\end{note}
\section{一般的说明}
估计程序对具有如下特征的动态模型是理想的：
\begin{itemize}
	\item 独立于模型。动态模型不遵循固定的结构形式(比如线性性)，独立于模型结构的估计程序最有用。
	\item 置信区间有解析形式。因为优化非线性动态模型以及自助所需复制的计算成本使得置信区间拥有解析形式是重要的。
	\item 误差项不要有太多假设。对动态模型误差项的独立同分布假设并不一直容易判断它是正确的，因此，不偏好该类假设的估计方法更佳。
	\item 多种数据类型均可使用。时序和截面应该都可以使用该方法。
\end{itemize}

本章我们提供模拟矩估计在动态模型中的应用的一个介绍。该法基本的思路就是定义数据合适的矩，通过改变不确定的参数，最小化这些矩和他们来自模型模拟的对应的矩的差距。矩是使用一些函数计算得到的统计量，这些函数以经验数据作为输入，并输出\textbf{一个}标量值。例如总体均值就是一个矩，以个体数据点为输入，以总体均值作为输出。任何函数都可以创造矩，当解析置信区间可行时，人们需要这些矩服从正态分布，这通常意味着每个矩都是来自统一潜在分布的多个独立观测值的函数的平均。实践上，这些观测值要么来自系统处于稳态时的时间序列(比如股价)，或者是生命中类似个体的类似观测点(比如某国所有五岁的个人)。典型地，这些个体具有不同的值，随机性具有重要作用(比如类似年龄个体的不同体重)。因此，存下如下问题特征时，MSM对动态建模有着最好的拟合：
\begin{itemize}
	\item 总体数据。MSM适合估计总体数据的模型估计。当诸如个人、公司、国家等不同数据单位都可行，对每个个体，一个或多个数据特征(如对于个人而言，包括体重、高度、年龄)是可行的。
	\item 随机过程的角色。模型包含驱动了模型的随机过程时，他们对模型行为的影响可以通过对参数的校准在数据中反映出来。(比如我们会尝试匹配由多个个体所观测到的方差。)
	\item 截面数据。MSM可以用于截面和时序数据。然而时序数据包含同一个体不同时间的多个数据点，截面数据包含同一时间多个个体的数据点。当数据是截面数据时，MSM也许是估计动态模型的唯一可行选择，因为它允许我们抽取隐藏在截面分布中的个体的历史轨迹。
	\item 置信区间。当寻找具有解析解的置信区间时，MSM是个合适的选择。
\end{itemize}

一个简单的关于个人身体重量的模型，身体重量由脂肪重量和非脂肪重量构成。缺乏时序数据，使用的是国家健康与营养调查数据库2005-2006，包含了5971个个体，这些个体基于种族(5个)、性别(2个)、年龄(11个)可以分为110个子组。在每个总体分组中，2个矩，体重的均值和方差被计算以匹配，这就会有220个矩用来匹配。

另一方面，该模型被5971个匹配了NHANES2006的人口学特征的个体所复制。初始的体重和脂肪比例是从另一个NHASNES1999-2000的样本中抽取的。注意到，NHANES的每轮调查都使用了不同其他轮的样板。因此，我们无法追踪时间上的同一个体，数据是截面的。模型进行了模拟，使得这个合成的总体从他们在2000年的最初年龄直到2006的最终年龄。在模拟的总体中不同子组的体重的均值和方差在2006年进行了计算，并和来自数据的220个矩进行了比较。利用每个矩的方差都倒数作为权重计算了误差平方的加权和，。。。该误差在Vensim内部优化引擎通过变化17个参数以最小化。估计的参数提供了最小的误差。因此，作者可以从包含不同年龄组的截面数据中估计一个包括个体成长机制都动态模型。

这个应用例子遵循了MSM的基本思想，但也不同于标准的MSM的程序 。首先，该应用中220个矩的数目高于许多典型的应用案例，一般而言，矩和待估参数的数目有着同样的数量级。第二，由于该设置的计算成本所限，在和数据进行比较前，每个矩仅模拟了一次，然而一般是使用不同的噪声种子模拟多次，然后提供对这个矩的估计。最终，该应用中并未报告置信区间。

\section{正式的定义}
假设经验数据在$ T $时刻，有$ N $个不同的观测值$ \{y_i\}_t,i\in\{1,\cdots,N\},t\in\{1,\cdots,T\} $。首先，需要计算经验矩。对于每个个体的总体组别，不同的矩条件可以被计算。这些矩函数放在一起形成一个$ p\times 1 $维的经验矩向量$ M_D $，该向量的第$ p $个典型元素可以记为$ m_p:m_p(\{y\}_t) $。

比如，我们可以说$ t $期$ N $个个体的一阶距条件是均值，$ M_D $中的对应元素为$ m_p(\{y\}_t)=\frac{1}{N}\sum_{i=1}^{N}\{y_i\}_t $。

接下来，需要确定模拟矩$ \tilde{M}_S $以匹配$ M_D $。一个完全设定的模型，即可以被一组参数所模拟，假设存在$ d $个元素的参数向量$ \theta $，系统动态的真实函数形式将产生$ \{y_i\}_t=g(\cdot) $，它可以被模型的输出$ \{\hat{y}_i\}_t=\hat{g}(\cdot) $所逼近。模型输出是已知参数$ C $、未知参数$ \theta $和随机冲击$ \xi $的函数。选择不同的$ \xi $将创造不同的$ \hat{g} $值，我们假设模型被正确设置，因此$ \hat{g}(\cdot) $是真实模型$g(\cdot) $的一致估计，即
\[ E(\hat{g}(C,\theta,\xi))=E(g(\cdot)) \]

这保证了如果我们使用真实的随机冲击流$ \xi $创造了足够多的样本输出，比如$ K  $个，模型输出的样本均值将创造与真实世界创造样本的一个合理逼近，即，
\[ m_p(\hat{g}(\cdot))=\frac{1}{K}\sum_{n=1}^Km_p(\hat{g}(C,\theta,\xi^{(n)})) \]

$ \tilde{M}_S $是模拟矩，它的每一个元素都由上式构造，以对应于经验矩$ M_D $。上式中的平均是非常重要的，这消除了样本偏差。一般相对小的$ K $足够保证精确的估计。然而，我们仍然需要足够大的重复以估计MSM的渐近性质。

一个识别模型的必要条件就是$ p\ge d $。未知参数可以如下估计，
\[ \hat{\theta}=\arg\min (\tilde{M}_S-M_D)'W (\tilde{M}_S-M_D)\]

很明显，要估计出$ d $个未知参数的必要(但非充分)条件是矩的数目$ p $要大于或等于待估参数的数目$ d $。大于的情况下成为过度识别，过度识别有一个检验。

MSM的核心思想就是样本矩和模拟矩之间的差异要最小。即通过选择$ \bm{\theta} $来最小化下式，
\[ \min_{\bm{\theta}}\quad (\bm{M}_s-\bm{M}_d)'\bm{W}(\bm{M}_s-\bm{M}_d) \]

其中，$ \bm{W} $是加权矩阵，后续我们会继续讨论加权矩阵的选择，本质上和计量上的广义矩中的加权矩阵是类似的。

\section{MSM的详细步骤}
\paragraph{第一步，选择矩条件。}一般模型输出的一阶和二阶矩是合适的候选者。最富有信息含量的矩应该：
\begin{itemize}
	\item 敏感于未知参数中的至少一个。	
	\item 有着相对小的方差，因为矩的方差太大，意味着噪声太多，因此是乏信息的。
\end{itemize}

此外除了单变量的一阶和二阶矩，也应该尝试变量间的矩，比如协方差，自相关，偏度等。

识别条件。选择正确的矩条件是模型识别和参数恢复最重要的一步。使用MSM的模型识别要求当且仅当结构参数等于真值时模型创造的矩条件匹配他们对应的经验部分。否则模型将创造出可疑的结果。而且，识别的充分条件就是结构参数与同维矩约束子集间的一对一映射。因为动态模型通常不会产生这样的闭式映射，为确保一个可识别的模型，人们应该选择一些敏感于结构参数的矩。检查这种敏感性的一个方法就是对每个被估参数做一个单维敏感性分析，以考察选择的矩条件是否随着参数的变化而变化。如果矩条件变化不大，或者它对参数变化的响应并不平滑和单调，那么这个矩也许并不是一个有信息含量的矩，而且会使得优化引擎在一个局部最优停下来，或者根本不会收敛到一个最优解。

\paragraph{第二步，加权矩阵}因为$ p $个矩条件，$ d $个未知参数，如果$ p=d $就是精确识别，$ p>d $就是过度识别，此时要使得样本矩和模拟矩之间的差异最小，OLS的思路又出来了。但在加总这个差异的时候，需要一个权重。类似于广义矩估计中的两步法，\textbf{首先}，选一个简单矩阵比如单位阵作为加权矩阵的初值。在计量经济学中通常是这么做的，但如果矩的量纲不同，使用单位阵还是有问题的。比如一个矩以千为单位，另一个矩以分数为单位，在优化中第一个矩会起主要作用，第二个矩会被忽略掉，这可能并不是我们想要的。因此，在这一步，我们需要把矩做一个尺度缩放。在这个思想指导下，就意味着在第一步我们不要用单位阵，而是用一个对角阵，这个对角阵的对角元素为$ 1/(M_d)^2 $来做我们第一步的加权矩阵$ \bm{W} $。

\textbf{然后，}来计算这个最优加权矩阵。利用上面讲的那个简单对角加权矩阵，模型中的未知参数$ \bm{\theta} $是可以估计出来的。这样，就可以通过重复模拟来得到模拟矩的方差协方差矩阵$ \hat{\bm{S}} $，
\begin{align*}
	\hat{\bm{S}}=\frac{1}{L_1}\sum_{l_1=1}^{L_1}\left[\bm{m}^{l_1}(\hat{g}(C,\hat{\bm{\theta}},\xi^{l_1}))-\frac{1}{L_2}\sum_{l_2=1}^{L_2}\bm{m}^{l_2}(\hat{g}(C,\hat{\bm{\theta}},\xi^{l_2}))\right]\cdot \\	
	\left[\bm{m}^{l_1}(\hat{g}(C,\hat{\bm{\theta}},\xi^{l_1}))-\frac{1}{L_2}\sum_{l_2=1}^{L_2}\bm{m}^{l_2}(\hat{g}(C,\hat{\bm{\theta}},\xi^{l_2}))\right]'
\end{align*}

其中，$ \bm{m}^{l_1} $表示利用系统第$ l_1 $次模拟出来的$ y_{it} $计算得到的$ p $维的模拟矩向量。
注意到没有$ \hat{\bm{\theta}} $的值，$ \hat{\bm{S}} $是计算不出来的，这也是为什么要用简单矩阵做第一步运算的道理。然后用$ \hat{\bm{S}}^{-1} $作为最优加权矩阵$ \bm{W}^* $。

\paragraph{第三部，初值和模拟。}对于非均衡模拟，我们想要从一个尽可能接近现实经验样本的初值开始模拟，这会减少由于初值不同带来的模拟矩的误差。在这样的模型中，我们仅能使用对应于经验样本时间的数据点。例如，我们仅能使用模拟模型中10岁大的数据与经验样本中10岁大的数据进行对比。对于稳态模型，所有的动态都在稳态，我可以在由模拟创造的数据的使用方面做的更有效率。本质上，无需运行模型$ N $次，而只要使用同一次模拟$ N $个后续的时间数据点即可。因为模型处于随机稳态，时间不同的点的差异体现了结果的稳态分布。如果这个过程被持续下去，我们应该忽略掉早期的观测值，
动态模型中的参数何以得到？我们这里简单介绍模拟矩估计。

\paragraph{第四步\quad 优化}因为是数值优化，所以初值对结果就显得很重要。推荐尝试不同的初值。另外要注意的时，尝试不同初值时，随机生成的样本不要有变化，这可以通过设置随机种子而做到。

\paragraph{第五步\quad 解释和推断}$ \bm{\theta} $的概率分布

\paragraph{一些建议}可以看到模拟矩估计对矩条件的选择有大的自由度，只要矩条件可行，研究者选择哪个矩条件是充分自由的。这当然也是一个风险。因此，提供一个矩对参数敏感性的图是一种保证矩富含信息的方法。另外每增加一个参数，估计的代价都是较高的。因此，一般建议从其他途径比如文献、回归等得到尽可能多的参数值，把剩下的参数估计再交给MSM。

\section{应用}
故事：200个公司10年记录的样本，公司有着类似的结构和参数，但有着不同的冲击。公司需要资源$ R $以得到产出$ O $。资源和产出的关系由函数$ f(R) $表达，
